From e6df788d92a6e16047771e4acfcb60e19b243e49 Mon Sep 17 00:00:00 2001
From: Kevin Schneider <schneider.kev@outlook.de>
Date: Mon, 21 Jun 2021 18:46:32 +0200
Subject: [PATCH] #5: Add function to create significance matrices from
 TMEAResults

---
 src/TMEA.Dash/TMEA.Dash.fsproj |  2 +-
 src/TMEA/Frames.fs             | 73 +++++++++++++++++++++++++++++-----
 src/TMEA/Playground.fsx        | 10 ++++-
 3 files changed, 71 insertions(+), 14 deletions(-)

diff --git a/src/TMEA.Dash/TMEA.Dash.fsproj b/src/TMEA.Dash/TMEA.Dash.fsproj
index 77734f1..bad22e9 100644
--- a/src/TMEA.Dash/TMEA.Dash.fsproj
+++ b/src/TMEA.Dash/TMEA.Dash.fsproj
@@ -26,7 +26,7 @@
   </ItemGroup>
 
   <ItemGroup>
-    <Content Include="TestFiles\OntologyMap.txt" />
+    <EmbeddedResource Remove="TestFiles\**" />
   </ItemGroup>
 
   <ItemGroup>
diff --git a/src/TMEA/Frames.fs b/src/TMEA/Frames.fs
index 86dc151..02f8eb9 100644
--- a/src/TMEA/Frames.fs
+++ b/src/TMEA/Frames.fs
@@ -4,13 +4,13 @@ module Frames =
 
     open Deedle
 
-    let private optDefFloat (f:float opt)=
+    let internal optDefFloat (f:float opt)=
         if f.HasValue then f.Value else -1.
 
-    let private optDefInt (f:int opt)=
+    let internal optDefInt (f:int opt)=
         if f.HasValue then f.Value else 0
 
-    let createTMEACharacterizationTable minBinSize (termNameTransformation: string -> string) (tmeaCharacterizations:TMEACharacterization []) : Frame<(string*(string*int)),string> =
+    let internal createTMEACharacterizationTable minBinSize (termNameTransformation: string -> string) (tmeaCharacterizations:TMEACharacterization []) : Frame<(string*(string*int)),string> =
         tmeaCharacterizations
         |> Array.mapi (fun i desc ->
             let negFrame: Frame<string,string> =
@@ -105,16 +105,67 @@ module Frames =
 
     type TMEAResult with
     
-        static member toTMEACharacterizationFrame(tmeaRes:TMEAResult) = 
-            tmeaRes.Characterizations |> createTMEACharacterizationTable 0 id
-        
-        static member toTMEACharacterizationFrame(minBinSize:int,tmeaRes:TMEAResult) = 
-            tmeaRes.Characterizations |> createTMEACharacterizationTable minBinSize id
+        static member toTMEACharacterizationFrame(?MinBinSize:int, ?TermNameTransformation:string->string) = 
+            
+            let minBinSize = defaultArg MinBinSize 0
+            let termNameTransformation = defaultArg TermNameTransformation id
+
+            fun (tmeaRes: TMEAResult) ->
+                tmeaRes.Characterizations |> createTMEACharacterizationTable minBinSize termNameTransformation
         
-        static member toTMEACharacterizationFrame (minBinSize:int,termNameTransformation:string->string,tmeaRes:TMEAResult) = 
-            tmeaRes.Characterizations |> createTMEACharacterizationTable minBinSize termNameTransformation
 
-        static member toSignificanceMatrixFrame : unit = raise (System.NotImplementedException())
+        static member toSignificanceMatrixFrame (?UseBenjaminiHochberg:bool, ?Threshold:float, ?TermNameTransformation:string->string) =
+            
+            let useBenjaminiHochberg = defaultArg UseBenjaminiHochberg false
+            let threshold = defaultArg Threshold 0.05
+            let termNameTransformation = defaultArg TermNameTransformation id
+
+            fun (tmeaRes:TMEAResult) ->
+                tmeaRes.Characterizations
+                |> Array.mapi (fun cI c ->
+                    let pos = 
+                        c.PositiveDescriptor
+                        |> Array.map (fun d ->
+                            d.OntologyTerm => d.PValue
+                        )
+            
+                    let neg = 
+                        c.PositiveDescriptor
+                        |> Array.map (fun d ->
+                            d.OntologyTerm => d.PValue
+                        )
+            
+                    let keys, pos' =
+                        if useBenjaminiHochberg then
+                            neg
+                            |> FSharp.Stats.Testing.MultipleTesting.benjaminiHochbergFDRBy id
+                            |> Array.ofList
+                            |> Array.unzip
+                        else
+                            neg
+                            |> Array.unzip
+            
+                    let keys, neg' =
+                        if useBenjaminiHochberg then
+                            neg
+                            |> FSharp.Stats.Testing.MultipleTesting.benjaminiHochbergFDRBy id
+                            |> Array.ofList
+                            |> Array.unzip
+                        else
+                            neg
+                            |> Array.unzip
+            
+            
+                    $"C_{cI}" => (
+                        Array.map2 (fun p n -> if p < n then p else n) pos' neg'
+                        |> Array.map (fun lowestPValue -> if lowestPValue < threshold then 1 else 0)
+                        |> Array.zip keys
+                        |> series
+                    )
+                )
+                |> frame
+                |> Frame.fillMissingWith 0
+                |> Frame.mapRowKeys (fun (name) -> (termNameTransformation name) => name)
 
         static member toConstraintsFrame : unit = raise (System.NotImplementedException())
 
diff --git a/src/TMEA/Playground.fsx b/src/TMEA/Playground.fsx
index ba9b449..fbeee16 100644
--- a/src/TMEA/Playground.fsx
+++ b/src/TMEA/Playground.fsx
@@ -68,7 +68,7 @@ open Plotly.NET
 
 let characterizationFrame = 
     testTmeaRes
-    |> TMEAResult.toTMEACharacterizationFrame
+    |> TMEAResult.toTMEACharacterizationFrame()
 
 characterizationFrame.Print()
 
@@ -79,7 +79,7 @@ testTmeaRes
 |> TMEAResult.plotConstraintPotentialTimecourses(OmitBaselineState=true,InvertConstraints=[|1;2;3|], ConstraintCutoff=3)
 
 testTmeaRes
-|> TMEAResult.plotPotentialHeatmap(ConstraintCutoff=3,InvertConstraints=[|2|])
+|> TMEAResult.plotConstraintPotentialHeatmap(ConstraintCutoff=3,InvertConstraints=[|2|])
 
 testTmeaRes
 |> TMEAResult.plotFreeEnergyLandscape()
@@ -94,6 +94,12 @@ testTmeaRes
 |> TMEAResult.plotFASWeightDistribution("signalling.light",[1;2;3],true,0.05)
 
 
+
+testTmeaRes
+|> TMEAResult.toSignificanceMatrixFrame()
+|> fun f -> f.Print()
+
+
 //readDataFrame 
 //    "TranscriptIdentifier" 
 //    "\t"