Merge pull request #12 from BioFSharp/developer

Import tested code from ProteomIQon

Author: caroott

RELEASE_NOTES.md

@@ -1,3 +1,10 @@
+#### 0.2.1 (Released 2026-08-30)
+
+* SparsePeakArray: new module upstreamed from ProteomIQon (sparse binned peak vectors over PeakArray with dot product)
+* SearchDB: add Sequence+GlobalMod ModSequence lookup upstreamed from ProteomIQon (prepareSelectModsequenceBySequenceAndGMod, getThreadSafePeptideLookUpFromFileBySequenceAndGMod)
+* ProteinInference: adopt the newer ProteomIQon behavior - PSMInput reads the ModelScore column (was PercolatorScore), assignTranscriptsToGenes takes a tryParseProteinID function instead of a regex string, createProteinModelInfoFromEntry reads the GFF3 ID attribute (was Name) and tolerates unknown strand characters as Forward (fixes a runtime match failure), isGene/isRNA generalized over the GFF line type parameter
+* FDRControl: add the PEP value machinery upstreamed from ProteomIQon (getLogisticRegressionFunction, createTargetDecoyHis, calculatePEPValues, logitTransformPepValues, initCalculateLin taking a trace callback instead of an NLog logger)
+
 #### 0.2.0 (Released 2026-08-29)
 
 * Upgrade to BioFSharp 2.0.0 (bundles the former BioFSharp.IO package) and migrate the codebase over the breaking API changes

src/BioFSharp.Mz/BioFSharp.Mz.fsproj

@@ -33,6 +33,7 @@
     <Compile Include="AssemblyInfo.fs" />
     <Compile Include="Peak.fs" />
     <Compile Include="PeakArray.fs" />
+    <Compile Include="SparsePeakArray.fs" />
     <Compile Include="PeakList.fs" />
     <Compile Include="Fragmentation.fs" />
     <Compile Include="Caching.fs" />

src/BioFSharp.Mz/FDRControl.fs

@@ -1,5 +1,6 @@
 namespace BioFSharp.Mz
 
+open System
 open FSharp.Stats
 open FSharpAux
 open FSharp.Stats.Fitting
@@ -247,3 +248,114 @@ module FDRControl =
         // takes a score from the dataset and assigns it a q value
         let interpolation = Interpolation.LinearSpline.predict linearSplineCoeff
         interpolation
+
+    /// for given data, creates a logistic regression model and returns a mapping function for this model
+    let getLogisticRegressionFunction (x:vector) (y:vector) epsilon =
+        let alpha =
+            match FSharp.Stats.Fitting.LogisticRegression.Univariable.estimateAlpha epsilon x y with
+            | Some a -> a
+            | None -> failwith "Could not find an alpha for logistic regression of fdr data"
+        let weight = FSharp.Stats.Fitting.LogisticRegression.Univariable.coefficient epsilon alpha x y
+        FSharp.Stats.Fitting.LogisticRegression.Univariable.predict weight
+
+    /// Creates a Histogram based on a given score of a target/decoy dataset. Each bin contains the information of the total count, the decoy count and the median score.
+    /// (Bin, Count, DecoyCount, Median Score)
+    let createTargetDecoyHis bandwidth (isDecoy: 'a -> bool) (decoyScoreF: 'a -> float) (targetScoreF: 'a -> float) (data: 'a[]) =
+        let halfBw = bandwidth / 2.0
+        let scoreDecoyInfo =
+            data
+            |> Array.map (fun x ->
+                if isDecoy x then
+                    {|Score = decoyScoreF x; Decoy = true|}
+                else
+                    {|Score = targetScoreF x; Decoy = false|}
+            )
+        scoreDecoyInfo
+        |> Array.groupBy (fun x ->
+            floor (x.Score / bandwidth))
+        |> Array.map (fun (k,values) ->
+            let count = (Array.length(values))
+            let decoyCount = (values |> Array.filter (fun x -> x.Decoy = true) |> Array.length)
+            let medianScore = values |> Array.map (fun x -> x.Score) |> Array.median
+            // first part of the tuple only needed for debugging
+            if k < 0. then
+                ((k  * bandwidth) + halfBw, count, decoyCount, medianScore)
+            else
+                ((k + 1.) * bandwidth - halfBw, count, decoyCount, medianScore)
+        )
+
+    /// Calculates the PEP value based on the ratio of Decoys to targets at a given score
+    let calculatePEPValues (totalCountF: 'a -> float) (decoyCountF: 'a -> float) (scoreF: 'a -> float) (dataFreq: 'a[]) =
+        dataFreq
+        |> Array.map (fun x ->
+            scoreF x,(decoyCountF x)/(totalCountF x)
+        )
+        |> Array.sortBy fst
+        |> Array.toList
+
+    /// Logit transforms pep values (log10)
+    let logitTransformPepValues score pepVal  =
+        Array.zip score pepVal
+        // 0 and 1 are + and - infinity
+        |> Array.filter (fun (y,x) -> x <> 0. && x <> 1.)
+        |> Array.map (fun (score,pep) ->
+            score,
+            log10 (pep/(1.-pep))
+        )
+        |> Array.unzip
+
+    /// Calculates monotonized PEP values for a target/decoy dataset based on the decoy/target ratio. Entries are binned with a given bandwidth as intital estiamtor based on the scores.
+    /// Returns a function which maps from score to PEP value based on a fit of a linear function using linear regression. The linear regression is performed on the logit transformed
+    /// pep values. The fit focuses on the pep values centered aound the middle of the score distribution
+    let initCalculateLin (trace: string -> unit) bandwidth (isDecoy: 'a -> bool) (decoyScoreF: 'a -> float) (targetScoreF: 'a -> float) (data: 'a[]) =
+        let lowerScore, upperScore =
+            let decoy =
+                data
+                |> Array.filter isDecoy
+                |> Array.map decoyScoreF
+                |> Array.filter (fun x -> x < 0.)
+                |> Array.median
+            let target =
+                data
+                |> Array.filter (isDecoy >> not)
+                |> Array.map targetScoreF
+                |> Array.filter (fun x -> x > 0.)
+                |> Array.median
+            decoy, target
+        trace (sprintf "Lower Score: %f; Upper Score: %f" lowerScore upperScore)
+        let filteredData =
+            data
+            |> Array.filter (fun entry ->
+                if isDecoy entry then
+                    let score = decoyScoreF entry
+                    score >= lowerScore && score <= upperScore
+                else
+                    let score = targetScoreF entry
+                    score >= lowerScore && score <= upperScore
+            )
+        trace (sprintf "Initial Bandwidth: %f" bandwidth)
+        let fittingFunction, score, pep =
+            let xPointRange =
+                let min = Math.Min((Array.minBy targetScoreF filteredData) |> targetScoreF, (Array.minBy decoyScoreF filteredData) |> decoyScoreF)
+                let max = Math.Max((Array.maxBy targetScoreF filteredData) |> targetScoreF, (Array.maxBy decoyScoreF filteredData) |> decoyScoreF)
+                max-min
+            let upperBW = Math.Min(10., xPointRange/10.)
+            [|bandwidth .. 0.1 .. upperBW|]
+            |> Array.choose (fun bw ->
+                let targetDecoyHis = createTargetDecoyHis bw (isDecoy: 'a -> bool) (decoyScoreF: 'a -> float) (targetScoreF: 'a -> float) (filteredData: 'a[])
+                let score',pep' =
+                    calculatePEPValues (fun (_,count,_,_) -> float count) (fun (_,_,decoyCount,_) -> float decoyCount) (fun (_,_,_,medianScore) -> medianScore) targetDecoyHis
+                    |> Array.ofList
+                    |> Array.unzip
+                let logitScore, logitPEPVal = logitTransformPepValues score' pep'
+                let coeff = Fitting.LinearRegression.OLS.Linear.Univariable.fit (vector logitScore) (vector logitPEPVal)
+                let fittingFunction' = (Fitting.LinearRegression.OLS.Linear.Univariable.predict coeff) >> (fun x -> 10.**(x)/(1.+10.**(x)))
+                let sos = FSharp.Stats.Fitting.GoodnessOfFit.calculateSumOfSquares fittingFunction' score' pep'
+                if coeff.[1] < 0. then
+                    Some (sos.Error/sos.Count, fittingFunction', score', pep', bw)
+                else
+                    None
+            )
+            |> Array.minBy (fun (error,_,_,_,_) -> error)
+            |> fun (error, fit,s,p,bw) -> trace (sprintf "Chosen Bandwidth: %f" bw); fit,s,p
+        fittingFunction

src/BioFSharp.Mz/ProteinInference.fs

@@ -73,7 +73,7 @@ module ProteinInference =
             PepSequenceID   : int
             [<FieldAttribute("StringSequence")>]
             Seq             :string
-            [<FieldAttribute("PercolatorScore")>]
+            [<FieldAttribute("ModelScore")>]
             Score : float
         }
 
@@ -165,13 +165,13 @@ module ProteinInference =
         String.filter (fun c -> System.Char.IsLower c |> not && c <> '[' && c <> ']') pepSeq
 
     /// Checks if GFF line describes gene
-    let isGene (item: GFFLine<seq<char>>) =
+    let isGene (item: GFFLine<'a>) =
         match item with
         | GFFEntryLine x -> x.Feature = "gene"
         | _ -> false
 
     /// Checks if GFF line describes rna
-    let isRNA (item: GFFLine<seq<char>>) =
+    let isRNA (item: GFFLine<'a>) =
         match item with
         | GFFEntryLine x -> if x.Feature = "mRNA" then Some x else None
         | _ -> None
@@ -183,7 +183,7 @@ module ProteinInference =
         let attributes = entry.Attributes
         /// Same as in FastA file
         let spliceVariantID =
-            match Map.tryFind "Name" attributes with
+            match Map.tryFind "ID" attributes with
             | Some res ->
                 res.Head
             | None ->
@@ -194,12 +194,13 @@ module ProteinInference =
             match entry.Strand with
             |'+' -> StrandDirection.Forward
             |'-' -> StrandDirection.Reverse
+            | _  -> StrandDirection.Forward
 
         createProteinModelInfo spliceVariantID chromosomeID direction locus i Seq.empty Seq.empty
 
     /// By reading GFF creates the protein models (relationships of proteins to each other) which basically means grouping the rnas over the gene loci
     /// TODO: Don't group over order but rather group over id
-    let assignTranscriptsToGenes regexPattern (gffLines: seq<GFFLine<seq<char>>>)  =
+    let assignTranscriptsToGenes tryParseProteinID (gffLines: seq<GFFLine<'a>>)  =
         gffLines
         // transcripts are grouped by the gene they originate from
         |> Seq.groupWhen isGene
@@ -212,13 +213,14 @@ module ProteinInference =
                 |> Seq.mapi (fun i element ->
                     // every transcript of gene gets its own number i and other info is collected from element and used for info of protein
                     let modelInfo = createProteinModelInfoFromEntry i locus element
-                    let r = System.Text.RegularExpressions.Regex.Match(modelInfo.Id,regexPattern)
+                    let r =  tryParseProteinID modelInfo.Id
                     // the gff3 id has to be matched with the sequence in the fasta file. therefore the regexpattern is used
-                    if r.Success then
-                        r.Value,
+                    match r with
+                    | Some v ->
+                        v,
                         modelInfo
-                    else
-                        failwithf "could not match gff3 entry id %s with regexpattern %s. Either gff3 file is corrupt or regexpattern is not correct" modelInfo.Id regexPattern
+                    | None ->
+                        failwithf "could not match gff3 entry id %s with regex pattern. Either gff3 file is corrupt or regexpattern is not correct" modelInfo.Id
                 )
 
             | _ -> Seq.empty

src/BioFSharp.Mz/SearchDB.fs

@@ -870,6 +870,20 @@ module SearchDB =
                 | true  -> Some (reader.GetDouble(0))
                 | false -> Option.None)
 
+            /// Prepares statement to select a ModSequence entry by Sequence and GlobalMod
+            let prepareSelectModsequenceBySequenceAndGMod (cn:SQLiteConnection) =
+                let querystring = "SELECT * FROM ModSequence WHERE Sequence=@sequence AND GlobalMod=@globalMod"
+                let cmd = new SQLiteCommand(querystring, cn)
+                cmd.Parameters.Add("@sequence", Data.DbType.String) |> ignore
+                cmd.Parameters.Add("@globalMod", Data.DbType.Int32) |> ignore
+                fun (sequence:string) (globalMod:int) ->
+                    cmd.Parameters.["@sequence"].Value  <- sequence
+                    cmd.Parameters.["@globalMod"].Value <- globalMod
+                    use reader = cmd.ExecuteReader()
+                    match reader.Read() with
+                        | true ->  (reader.GetInt32(0), reader.GetInt32(1),reader.GetDouble(2), reader.GetInt64(3), reader.GetString(4), reader.GetInt32(5))
+                        | false -> -1,-1,nan,-1L,"",-1
+
             /// Prepares statement to select a Protein Accession entry by ID
             let prepareSelectProteinAccessionByID (cn:SQLiteConnection) (tr) =
                 let querystring = "SELECT Accession FROM Protein WHERE ID=@id "
@@ -1671,7 +1685,15 @@ module SearchDB =
                 selectModsequenceByMassRange lowerMass' upperMass'
                 |> List.map (createLookUpResultBy parseAAString)
         )   
-   
+
+    /// Returns a LookUpResult
+    let getThreadSafePeptideLookUpFromFileBySequenceAndGMod (cn:SQLiteConnection) sdbParams =
+        let parseAAString = initOfModAminoAcidString sdbParams.IsotopicMod (sdbParams.FixedMods@sdbParams.VariableMods)
+        let selectModsequenceByID = Db.SQLiteQuery.prepareSelectModsequenceBySequenceAndGMod cn
+        (fun sequence globalMod ->
+            selectModsequenceByID sequence globalMod
+            |> (createLookUpResultBy parseAAString))
+
     let copyDBIntoMemory (cn:SQLiteConnection) = 
         //cn.Open()
         let inMemoryDB = new SQLiteConnection("Data Source=:memory:;cache=shared;Version=3")

src/BioFSharp.Mz/SparsePeakArray.fs

@@ -0,0 +1,46 @@
+namespace BioFSharp.Mz
+
+module SparsePeakArray =
+
+    ///
+    type SparsePeakArray = {
+        Data            : System.Collections.Generic.IDictionary<int,float>
+        MzToBinIdx      : float -> int
+        BinIdxToMz      : int -> float
+        }
+
+    ///
+    let dot (x:SparsePeakArray) (y:SparsePeakArray) =
+        x.Data
+        |> Seq.fold (fun (acc:float) xi ->
+            let present,yi = y.Data.TryGetValue xi.Key
+            if present then acc + (yi * xi.Value) else acc
+            ) 0.
+
+    ///
+    let initMzToBinIdx width offset x = int ((x / width) + offset)
+
+    let initBinIdxToMz width offset x =
+        ((float x) - offset) * width
+
+
+    ///
+    let peaksToNearestBinVector binWidth offset (minMassBoarder:float) (maxMassBoarder:float) (pkarr:PeakArray<_>) =
+        let mzToBinIdx = initMzToBinIdx binWidth offset
+        let binIdxToMz = initBinIdxToMz binWidth offset
+        let keyValues =
+            pkarr
+            |> Array.choose (fun p ->
+                if p.Mz < maxMassBoarder && p.Mz > minMassBoarder then
+                    let index = (mzToBinIdx p.Mz)
+                    Some (index, p.Intensity)
+                else
+                    None
+                )
+            |> Array.groupBy fst
+            |> Array.map (fun (idx,data) -> idx, data |> Array.sumBy snd)
+        {
+            Data        = keyValues |> dict
+            MzToBinIdx  = mzToBinIdx
+            BinIdxToMz  = binIdxToMz
+        }

tests/BioFSharp.Mz.Tests/BioFSharp.Mz.Tests.fsproj

@@ -8,6 +8,7 @@
 
   <ItemGroup>
     <Compile Include="PeakTests.fs" />
+    <Compile Include="SparsePeakArrayTests.fs" />
     <Compile Include="CachingTests.fs" />
     <Compile Include="StatsExtensionTests.fs" />
     <Compile Include="SignalDetectionTests.fs" />