Port existing modules

Author: kMutagene

src/BioFSharp.Stats/BioFSharp.Stats.fsproj

@@ -6,11 +6,15 @@
   </PropertyGroup>
 
   <ItemGroup>
-    <Compile Include="Library.fs" />
+    <Compile Include="PhylTree.fs" />
+    <Compile Include="OntologyEnrichment.fs" />
+    <Compile Include="RNASeq.fs" />
+    <Compile Include="SAEPT.fs" />
   </ItemGroup>
 
   <ItemGroup>
     <PackageReference Include="BioFSharp" Version="2.0.0-preview.3" />
+    <PackageReference Include="FSharp.Stats" Version="0.4.11" />
   </ItemGroup>
 
 </Project>

src/BioFSharp.Stats/Library.fs

@@ -0,0 +1,161 @@
+namespace BioFSharp.Stats
+
+open System
+open FSharpAux
+
+module OntologyEnrichment =
+
+    /// Represents an item in an ontology set
+    type OntologyItem<'a> = {
+        Id               : string
+        OntologyTerm     : string
+        GroupIndex       : int
+        Item             : 'a
+        }
+
+
+    /// Creates an item in an ontology set
+    let createOntologyItem id ontologyTerm groupIndex item =
+        {Id = id; OntologyTerm = ontologyTerm; GroupIndex = groupIndex; Item = item}
+
+    /// Represents a gene set enrichment result 
+    type GseaResult<'a> = {
+        ///Ontology term e.g. MapMan term, GO term ...
+        OntologyTerm     : string
+        ///Sequence of single items associated to the ontology term 
+        ItemsInBin       : seq<OntologyItem<'a>>
+        ///Number of significantly altered items in 'OntologyTerm' bin
+        NumberOfDEsInBin : int
+        ///Number of items in 'OntologyTerm' bin
+        NumberInBin      : int
+        ///Number of significantly altered items within the total data set
+        TotalNumberOfDE  : int
+        ///Number of all items (expanded)
+        TotalUniverse    : int
+        PValue           : float
+        }
+
+    /// Creates a gene set enrichment result 
+    let createGseaResult ontologyTerm desInBin numberOfDEsInBin numberInBin totalNumberOfDE totalUnivers pValue = 
+        {OntologyTerm = ontologyTerm;ItemsInBin = desInBin; NumberOfDEsInBin = numberOfDEsInBin; 
+            NumberInBin = numberInBin; TotalNumberOfDE = totalNumberOfDE; TotalUniverse = totalUnivers; PValue = pValue}
+
+    ///Splits an OntologyEntry with seperator concatenated TermIds
+    let splitMultipleAnnotationsBy (separator:char) (item:OntologyItem<'A>) =
+        let annotations = item.OntologyTerm.Split(separator)
+        annotations
+        |> Seq.map (fun ot -> {item with OntologyTerm = ot})
+
+    /// Splits MapMan OntologyEntries with seperator concatenated TermIds
+    /// Attention: Also parses string to int to get rid of 0 - terms
+    let splitMapManOntologyItemsBy (separator:char) (data:seq<OntologyItem<'a>>) =
+        let splitTerm (termId:string) (separator:char) =
+            termId.Split(separator) 
+            |> Array.map (fun sTerm -> 
+                let splited = sTerm.Split('.')
+                let toInt = splited |> Seq.map (fun v -> Int32.Parse(v).ToString())                                                                            
+                toInt  |> String.concat "." 
+                         )
+        data
+        |> Seq.collect (fun oi -> 
+            splitTerm oi.OntologyTerm separator
+            |> Seq.map (fun sTerm -> createOntologyItem oi.Id sTerm oi.GroupIndex oi.Item)
+                        )
+
+
+    /// Extends leaf OntologyEntries to their full tree
+    let expandOntologyTree (data:seq<OntologyItem<'a>>) =
+        data
+        |> Seq.collect (fun oi -> 
+            let expandenTermIds = oi.OntologyTerm.Split('.') |> Array.scanReduce (fun acc elem -> acc + "." + elem)
+            expandenTermIds |> Seq.map (fun sTerm -> createOntologyItem oi.Id sTerm oi.GroupIndex oi.Item) 
+                       )
+
+
+    // ###########################################################################################################
+    // the hypergeometric distribution is a discrete probability distribution that describes the probability of 
+    //   k successes in
+    //   n draws from a finite 
+    //   x population of size containing
+    //   m successes without replacement (successes states)
+    /// Calculates p value based on hypergeometric distribution (pValue <= k)
+    let CalcHyperGeoPvalue numberOfDEsInBin numberInBin totalUnivers totalNumberOfDE (splitPvalueThreshold:int) =
+        if (numberOfDEsInBin > 1) then
+            let hp = FSharp.Stats.Distributions.Discrete.Hypergeometric.Init totalUnivers totalNumberOfDE numberInBin            
+            if numberInBin > splitPvalueThreshold then                                
+                // Calculate normal pValue
+                1. -  hp.CDF (float (numberOfDEsInBin - 1)) 
+            else
+                // Calculate split pValue
+                0.5 * ((1. -  hp.CDF(float(numberOfDEsInBin - 1)) ) + ( (1. -  hp.CDF(float(numberOfDEsInBin))) ) )
+        else
+                nan
+        
+       
+    // #######################################################    
+    //  functional term enrichment is calculated according to following publication
+    //  http://bioinformatics.oxfordjournals.org/cgi/content/abstract/23/4/401
+    //  also includes mid-pValues
+    /// Calculates functional term enrichment
+    let CalcSimpleOverEnrichment (deGroupIndex:int) (splitPvalueThreshold:option<int>) (data:seq<OntologyItem<'a>>) =
+        let _splitPvalueThreshold   = defaultArg splitPvalueThreshold 5           
+        
+        let totalUnivers    = data |> Seq.length
+        let totalNumberOfDE = data |> Seq.filter (fun oi -> oi.GroupIndex = deGroupIndex) |> Seq.length
+        
+        // returns (DE count, all count)
+        let countDE (subSet:seq<OntologyItem<'a>>) =             
+            let countMap = 
+                subSet 
+                |> Seq.countBy (fun oi -> if oi.GroupIndex = deGroupIndex then true else false )
+                |> Map.ofSeq
+            (countMap.TryFindDefault 0 true,(countMap.TryFindDefault 0 true) + (countMap.TryFindDefault 0 false))
+        
+        data
+        |> Seq.groupBy ( fun oi -> oi.OntologyTerm)
+        |> Seq.map (fun (oTerm,values) -> 
+            let numberOfDEsInBin,numberInBin = countDE values
+            let pValue = CalcHyperGeoPvalue numberOfDEsInBin numberInBin totalUnivers totalNumberOfDE _splitPvalueThreshold
+            createGseaResult oTerm values numberOfDEsInBin numberInBin totalNumberOfDE totalUnivers pValue)
+
+
+    // #######################################################    
+    //  functional term enrichment is calculated according to following publication
+    //  http://bioinformatics.oxfordjournals.org/cgi/content/abstract/23/4/401
+    //  also includes mid-pValues
+    /// Calculates functional term enrichment
+    let CalcOverEnrichment (deGroupIndex:int) (splitPvalueThreshold:option<int>) (minNumberInTerm:option<int>) (data:seq<OntologyItem<'a>>) =
+        let _splitPvalueThreshold   = defaultArg splitPvalueThreshold 5
+        let _minNumberInTerm        = defaultArg minNumberInTerm 2
+        
+        // Distinct by term and gene name
+        // Has to be done by an ouside function
+        //let distinctData    = data |> Seq.distinctBy (fun o -> o.displayID)                
+        let gData           = data |> Seq.groupBy ( fun o -> o.OntologyTerm)
+        // reduce to terms at least annotated with 2 items
+        let fData = gData |> Seq.filter ( fun (key:string,values:seq<OntologyItem<'a>>) -> Seq.length(values) >= _minNumberInTerm)
+        let groupCount = fData |> Seq.collect (fun (key:string,values:seq<OntologyItem<'a>>) -> values ) |> Seq.countBy (fun o -> o.GroupIndex)
+        
+        let totalUnivers    = groupCount |> Seq.fold (fun  (acc:int) (index:int,count:int) -> acc + count) 0
+        let totalNumberOfDE = 
+            let tmp = groupCount |> Seq.tryFind (fun (key,v) -> key = deGroupIndex)
+            if tmp.IsNone then 
+                raise (System.ArgumentException("DE group index does not exists in ontology entry"))
+            else
+                snd(tmp.Value)
+        
+        // returns (DE count, all count)
+        let countDE (subSet:seq<OntologyItem<'a>>) =             
+            let countMap = 
+                subSet 
+                |> Seq.countBy (fun (oi) -> oi.GroupIndex = deGroupIndex)
+                |> Map.ofSeq
+            (countMap.TryFindDefault 0 true,(countMap.TryFindDefault 0 true) + (countMap.TryFindDefault 0 false))
+        
+        fData
+        |> Seq.map (fun (oTerm,values) -> 
+            let numberOfDEsInBin,numberInBin = countDE values
+            let pValue = CalcHyperGeoPvalue numberOfDEsInBin numberInBin totalUnivers totalNumberOfDE _splitPvalueThreshold
+            createGseaResult oTerm values numberOfDEsInBin numberInBin totalNumberOfDE totalUnivers pValue)
+
+

src/BioFSharp.Stats/OntologyEnrichment.fs

@@ -0,0 +1,48 @@
+namespace BioFSharp.Stats
+
+open BioFSharp
+open FSharp.Stats.ML.Unsupervised
+
+module PhylTree =
+
+    /// <summary>
+    /// converts the input hierarchical clustering to a phylogenetig tree and conserves the distance insformation.
+    /// In contrast to the clustering result, the distance value of a Branch represents the distance to its Parent, 
+    /// not the distance that all children have to this Branch.
+    /// </summary>
+    let ofHierarchicalCluster (branchTag:'T) (distanceConverter: float -> 'Distance) (hCluster:HierarchicalClustering.Cluster<'T>) : PhylogeneticTree<'T * 'Distance>=
+        let rec loop distance (c: HierarchicalClustering.Cluster<'T>) =
+            match c with
+            | HierarchicalClustering.Cluster.Node (cIndex, distance, lCount, left, right) ->    
+                PhylogeneticTree.Branch ((branchTag, distanceConverter distance), [loop distance left; loop distance right])
+            | HierarchicalClustering.Cluster.Leaf (id, lCount, tag) -> PhylogeneticTree.Leaf((tag, distanceConverter distance))
+        loop 0. hCluster 
+
+    /// <summary>
+    /// Performs hierarchical clustering of the input TaggedSequences using the provided distance function and linker. Returns the result as a Phylogenetic tree.</summary>
+    /// </summary>
+    /// <parameter name="branchTag">a tag to give the infered common ancestor branches (these are not tagged in contrast to the input sequence.)</parameter>
+    /// <parameter name="distanceConverter">a converter function for the distance between nodes of the tree. Usually, a conversion to a string makes sense for downstream conversion to Newick format</parameter>
+    /// <parameter name="distanceFunction">a function that determines the distance between two sequences e.g. evolutionary distance based on a substitution model</parameter>
+    /// <parameter name="linker">the linker function to join clusters with</parameter>
+    /// <parameter name="sequences">the input TaggedSequences</parameter>
+    let ofTaggedSequencesWithLinker (branchTag:'T) (distanceConverter: float -> 'Distance) (distanceFunction: seq<'S> -> seq<'S> -> float) linker (sequences: seq<TaggedSequence<'T,'S>>) =
+        sequences
+        |> HierarchicalClustering.generate
+            (fun a b  -> distanceFunction a.Sequence b.Sequence)
+            linker
+        |> ofHierarchicalCluster (TaggedSequence.create branchTag Seq.empty) distanceConverter
+
+
+    /// <summary>
+    /// Performs hierarchical clustering of the input TaggedSequences using the provided distance function. Returns the result as a Phylogenetic tree.</summary>
+    /// </summary>
+    /// <parameter name="distanceFunction">a function that determines the distance between two sequences e.g. evolutionary distance based on a substitution model</parameter>
+    /// <parameter name="sequences">the input TaggedSequences</parameter>
+    let ofTaggedBioSequences (distanceFunction: seq<#IBioItem> -> seq<#IBioItem> -> float) (sequences: seq<TaggedSequence<string,#IBioItem>>) : PhylogeneticTree<TaggedSequence<string,#IBioItem>*float>=
+        sequences
+        |> ofTaggedSequencesWithLinker
+            "Ancestor"
+            id
+            distanceFunction
+            HierarchicalClustering.Linker.upgmaLwLinker

src/BioFSharp.Stats/PhylTree.fs

@@ -0,0 +1,87 @@
+namespace BioFSharp.Stats
+
+/// Contains types and functions needed for RNA-Seq normalization
+module RNASeq = 
+
+    /// Input type for RNA-Seq normalization
+    type RNASeqInput = { 
+        GeneID : string
+        GeneLength : float
+        GeneCount : float
+    } with
+        static member Create(
+            id: string,
+            gl: float,
+            gc: float
+        ) = { GeneID=id; GeneLength=gl; GeneCount=gc}
+
+    type NormalizationMethod = 
+        | RPKM
+        | TPM
+
+    /// Type with GeneID, normalized data and method of normalization
+    type NormalizedCounts = { 
+        GeneID : string
+        NormalizedCount : float
+        NormalizationMethod: NormalizationMethod
+    } with 
+        static member Create(
+            id: string,
+            nc: float,
+            nm: NormalizationMethod
+        ) = { GeneID=id; NormalizedCount=nc; NormalizationMethod=nm}
+
+    /// calculates Reads Per Million
+    let private calcRPM sumOfAllReadsPerMil counts =
+        (counts |> float) / sumOfAllReadsPerMil
+
+    /// calculates RPKM
+    let private calcRPKM geneLength rpm = 
+        (float rpm) / ((float geneLength) / 1000.)
+
+    ///Performs RPKM normalization
+    let private rpkmsOf (geneIDs:seq<string>) (length:seq<float>) (counts:seq<float>) =
+
+        let sumOfAllReads = counts |> Seq.sum
+        let sumOfAllReadsPerMil = sumOfAllReads / 1000000. 
+        let rpms = Seq.map (fun counts -> calcRPM sumOfAllReadsPerMil counts) counts
+
+        let rpkms =
+            Seq.zip length rpms
+            |> Seq.map (fun (length, rpm) -> calcRPKM length rpm)
+
+        let rpkmResult =
+            Seq.map2 (fun ids counts -> NormalizedCounts.Create(ids, counts, RPKM)) geneIDs rpkms
+
+        rpkmResult
+
+    /// Returns RPKM normalized data
+    let rpkms (idLengthAndCounts:seq<RNASeqInput>) =
+        rpkmsOf 
+            (idLengthAndCounts |> Seq.map (fun x -> x.GeneID)) 
+            (idLengthAndCounts |> Seq.map (fun x -> x.GeneLength)) 
+            (idLengthAndCounts |> Seq.map (fun x -> x.GeneCount)) 
+    
+    /// Performs TPM normalization
+    let private tpmsOf (idLengthAndCounts:seq<RNASeqInput>) =
+        let rpk = 
+            idLengthAndCounts
+            |> Seq.map (fun idLengthAndCounts -> idLengthAndCounts.GeneCount/idLengthAndCounts.GeneLength/1000.)
+        
+        let sumOfAllReads = rpk |> Seq.sum
+        let sumOfAllReadsPerMil = sumOfAllReads / 1000000.
+        let tpms = rpk |> Seq.map (fun rpks -> rpks/sumOfAllReadsPerMil)
+        
+        let geneID =
+            idLengthAndCounts
+            |> Seq.map (fun idLengthAndCounts -> idLengthAndCounts.GeneID)
+        
+        let tpmResult =
+            Seq.map2 (fun ids counts -> NormalizedCounts.Create(ids, counts, TPM)) geneID tpms
+        
+        tpmResult
+
+    /// Returns TPM normalized data
+    let tpms (idLengthAndCounts:seq<RNASeqInput>) =
+        tpmsOf idLengthAndCounts 
+