Chisq estimation for /match

2 files changed, 92 insertions, 13 deletions

diff --git a/bbrc-sample b/bbrc-sample
-Subproject 6ddfc2dc414f1e64ac16286c0cee5a4b0022d2e
+Subproject 0d1d349ac33ae2fcc1bbdf31617ed9132c7527c
diff --git a/fminer.rb b/fminer.rb
index d8da725..9942cfa 100644
--- a/fminer.rb
+++ b/fminer.rb
@@ -119,30 +119,109 @@ get "/fminer/last/?" do
   end
 end
 
-# Creates same features for dataset <dataset_uri> that have been created
-# with fminer in dataset <feature_dataset_uri>
-# accept params[:nr_hits] as used in other fminer methods
+# Matches features of a a feature dataset onto instances of another dataset. 
+# The latter is referred to as 'training dataset', since p-values are computed,
+# if user passes a prediction feature, or if the training dataset has only one feature.
+# The result does not contain the prediction feature.
+# @param [String] dataset_uri URI of the dataset 
+# @param [String] feature_dataset_uri URI of the feature dataset (i.e. dependent variable)
+# @param [optional] parameters Accepted parameters are
+# - prediction_feature URI of prediction feature to calculate p-values for
+# @return [text/uri-list] Task URI
 post '/fminer/:method/match?' do
   raise OpenTox::BadRequestError.new "feature_dataset_uri not given" unless params[:feature_dataset_uri]
   raise OpenTox::BadRequestError.new "dataset_uri not given" unless params[:dataset_uri]
+
+  training_dataset = OpenTox::Dataset.find "#{params[:dataset_uri]}"
+  unless params[:prediction_feature] # try to read prediction_feature from dataset
+    prediction_feature = OpenTox::Feature.find(training_dataset.features.keys.first) if training_dataset.features.size == 1
+  end
+  prediction_feature = OpenTox::Feature.find(params[:prediction_feature]) if params[:prediction_feature]
+
   task = OpenTox::Task.create("Matching features", url_for('/fminer/match',:full)) do |task|
+
+    # get endpoint statistics
+    if prediction_feature
+      db_class_sizes = Array.new # for effect calculation
+      all_activities = Hash.new # for effect calculation, indexed by id, starting from 1 (not 0)
+      id = 1
+      training_dataset.compounds.each do |compound|
+        entry=training_dataset.data_entries[compound]
+        entry.each do |feature,values|
+          if feature == prediction_feature.uri
+            values.each { |val|
+              if val.nil? 
+                LOGGER.warn "No #{feature} activity for #{compound.to_s}."
+              else
+                if prediction_feature.feature_type == "classification"
+                  activity= training_dataset.value_map(prediction_feature.uri).invert[val].to_i # activities are mapped to 1..n
+                  db_class_sizes[activity-1].nil? ? db_class_sizes[activity-1]=1 : db_class_sizes[activity-1]+=1 # AM effect
+                elsif prediction_feature.feature_type == "regression"
+                  activity= val.to_f 
+                end
+                begin
+                  all_activities[id]=activity # DV: insert global information
+                  id += 1
+                rescue Exception => e
+                  LOGGER.warn "Could not add " + smiles + "\t" + val.to_s + " to fminer"
+                  LOGGER.warn e.backtrace
+                end
+              end
+            }
+          end
+        end
+      end
+    end
+
+    # Intialize result by adding compounds
     f_dataset = OpenTox::Dataset.find params[:feature_dataset_uri],@subjectid
     c_dataset = OpenTox::Dataset.find params[:dataset_uri],@subjectid
     res_dataset = OpenTox::Dataset.create CONFIG[:services]["dataset"],@subjectid
-    f_dataset.features.each do |f,m|
-      res_dataset.add_feature(f,m)
-    end
     c_dataset.compounds.each do |c|
       res_dataset.add_compound(c)
     end
+
+    # Run matching, put data entries in result. Features are recreated.
     smi = [nil]; smi += c_dataset.compounds.collect { |c| OpenTox::Compound.new(c).to_smiles }
     smarts = f_dataset.features.collect { |f,m| m[OT.smarts] }
     params[:nr_hits] == "true" ? hit_count=true: hit_count=false
-    matches, counts = LU.new.match_rb(smi, smarts, hit_count)
+    matches, counts = LU.new.match_rb(smi, smarts, hit_count) if smarts.size>0
+
     f_dataset.features.each do |f,m|
       if (matches[m[OT.smarts]] && matches[m[OT.smarts]].size>0)
+
+        feature_uri = File.join res_dataset.uri,"feature","match", res_dataset.features.size.to_s
+        metadata = {
+          RDF.type => [OT.Feature, OT.Substructure],
+          OT.hasSource => f_dataset.uri,
+          OT.smarts => m[OT.smarts],
+          OT.parameters => [
+            { DC.title => "dataset_uri", OT.paramValue => params[:dataset_uri] }
+          ]
+        }
+
+        if (prediction_feature) 
+          feat_hash = Hash[*(all_activities.select { |k,v| matches[m[OT.smarts]].include?(k) }.flatten)]
+          if prediction_feature.feature_type == "regression"
+            p_value = @@last.KSTest(all_activities.values, feat_hash.values).to_f # AM LAST: use internal function for test
+            effect = (p_value > 0) ? "activating" : "deactivating"
+          else
+            p_value = @@last.ChisqTest(all_activities.values, feat_hash.values).to_f
+            g=Array.new # g is filled in *a*scending activity
+            training_dataset.value_map(prediction_feature.uri).each { |y,act| g[y-1]=Array.new }
+            feat_hash.each  { |x,y|   g[y-1].push(x)   }
+            max = OpenTox::Algorithm.effect(g, db_class_sizes) # db_class_sizes is filled in *a*scending activity
+            effect = max+1
+          end
+          metadata[OT.effect] = effect
+          metadata[OT.pValue] = p_value.abs
+          metadata[OT.parameters] << { DC.title => "prediction_feature", OT.paramValue => prediction_feature.uri }
+        end
+        
+        res_dataset.add_feature feature_uri, metadata
+
         matches[m[OT.smarts]].each_with_index {|id,idx| 
-          res_dataset.add(c_dataset.compounds[id-1],f,counts[m[OT.smarts]][idx])
+          res_dataset.add(c_dataset.compounds[id-1],feature_uri,counts[m[OT.smarts]][idx])
         }
       end
     end
@@ -225,9 +304,9 @@ post '/fminer/bbrc/?' do
         p_value = f[1]
 
         if (!@@bbrc.GetRegression)
-          id_arrs = f[2..-1].flatten
-          max = OpenTox::Algorithm.effect(f[2..-1], fminer.db_class_sizes)
-          effect = f[2..-1].size-max
+          id_arrs = f[2..-1].flatten                                       # f[2..-1] is filled in *de*scending order,
+          max = OpenTox::Algorithm.effect(f[2..-1], fminer.db_class_sizes) # db_class_size is filled in *a*scending order,
+          effect = f[2..-1].size-max                                       # thus need to turn around effect
         else #regression part
           id_arrs = f[2]
           # DV: effect calculation
@@ -344,7 +423,7 @@ post '/fminer/bbrc/sample/?' do
 
   # method
   unless params[:method]
-    method="mean"
+    method="mle"
     LOGGER.debug "Set method to default value #{method}"
   else
     raise OpenTox::BadRequestError.new "method is neither 'mle' nor 'mean'" unless (params[:method] == "mle" or params[:method] == "mean")
@@ -436,7 +515,7 @@ post '/fminer/bbrc/sample/?' do
       @value_map.each { |y,act| g[y-1]=Array.new }
       feat_hash.each  { |x,y|   g[y-1].push(x)   }
       max = OpenTox::Algorithm.effect(g, fminer.db_class_sizes)
-      effect = g.size-max
+      effect = max + 1
       feature_uri = File.join feature_dataset.uri,"feature","bbrc", features.size.to_s
       unless features.include? smarts
         features << smarts