diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 82e0051..f63a9cd 100644 --- a/RELEASE_NOTES.md +++ b/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 diff --git a/src/BioFSharp.Mz/BioFSharp.Mz.fsproj b/src/BioFSharp.Mz/BioFSharp.Mz.fsproj index b04d3c6..e011efb 100644 --- a/src/BioFSharp.Mz/BioFSharp.Mz.fsproj +++ b/src/BioFSharp.Mz/BioFSharp.Mz.fsproj @@ -33,6 +33,7 @@ + diff --git a/src/BioFSharp.Mz/FDRControl.fs b/src/BioFSharp.Mz/FDRControl.fs index 6532bd0..95cfaf7 100644 --- a/src/BioFSharp.Mz/FDRControl.fs +++ b/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 diff --git a/src/BioFSharp.Mz/ProteinInference.fs b/src/BioFSharp.Mz/ProteinInference.fs index 9c8a501..41f6746 100644 --- a/src/BioFSharp.Mz/ProteinInference.fs +++ b/src/BioFSharp.Mz/ProteinInference.fs @@ -73,7 +73,7 @@ module ProteinInference = PepSequenceID : int [] Seq :string - [] + [] 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>) = + let isGene (item: GFFLine<'a>) = match item with | GFFEntryLine x -> x.Feature = "gene" | _ -> false /// Checks if GFF line describes rna - let isRNA (item: GFFLine>) = + 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>>) = + let assignTranscriptsToGenes tryParseProteinID (gffLines: seq>) = 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 diff --git a/src/BioFSharp.Mz/SearchDB.fs b/src/BioFSharp.Mz/SearchDB.fs index 88deadb..8aa1e0c 100644 --- a/src/BioFSharp.Mz/SearchDB.fs +++ b/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") diff --git a/src/BioFSharp.Mz/SparsePeakArray.fs b/src/BioFSharp.Mz/SparsePeakArray.fs new file mode 100644 index 0000000..1dd382a --- /dev/null +++ b/src/BioFSharp.Mz/SparsePeakArray.fs @@ -0,0 +1,46 @@ +namespace BioFSharp.Mz + +module SparsePeakArray = + + /// + type SparsePeakArray = { + Data : System.Collections.Generic.IDictionary + 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 + } diff --git a/tests/BioFSharp.Mz.Tests/BioFSharp.Mz.Tests.fsproj b/tests/BioFSharp.Mz.Tests/BioFSharp.Mz.Tests.fsproj index 9d3a7e5..3366c9d 100644 --- a/tests/BioFSharp.Mz.Tests/BioFSharp.Mz.Tests.fsproj +++ b/tests/BioFSharp.Mz.Tests/BioFSharp.Mz.Tests.fsproj @@ -8,6 +8,7 @@ + diff --git a/tests/BioFSharp.Mz.Tests/FDRControlTests.fs b/tests/BioFSharp.Mz.Tests/FDRControlTests.fs index 9b00c98..8c83ac0 100644 --- a/tests/BioFSharp.Mz.Tests/FDRControlTests.fs +++ b/tests/BioFSharp.Mz.Tests/FDRControlTests.fs @@ -2,6 +2,7 @@ module FDRControlTests open System open Expecto +open FSharp.Stats open BioFSharp.Mz let private expectWithin tolerance actual expected message = @@ -347,6 +348,93 @@ let tests = // domain reading of a perfectly separated target/decoy score distribution. ] + testList "PEP" [ + testCase "createTargetDecoyHis bins scores with half-bandwidth-centered labels and per-bin decoy counts" <| fun _ -> + let data = [| (0.2, false); (0.7, true); (0.4, false); (1.3, false) |] + let bins = + FDRControl.createTargetDecoyHis + 1.0 + snd + (fun (score, _) -> score) + (fun (score, _) -> score) + data + |> Array.sortBy (fun (bin, _, _, _) -> bin) + Expect.equal bins.Length 2 "the data is split into two score bins" + let bin0, count0, decoyCount0, median0 = bins.[0] + expectFloatClose bin0 0.5 "the non-negative bin label is centered by half the bandwidth" + Expect.equal count0 3 "the first bin contains three entries" + Expect.equal decoyCount0 1 "the first bin contains one decoy" + expectFloatClose median0 0.4 "the first bin median score" + let bin1, count1, decoyCount1, median1 = bins.[1] + expectFloatClose bin1 1.5 "the second non-negative bin label is centered by half the bandwidth" + Expect.equal count1 1 "the second bin contains one entry" + Expect.equal decoyCount1 0 "the second bin contains no decoys" + expectFloatClose median1 1.3 "the second bin median score" + let negativeBin, _, _, _ = + FDRControl.createTargetDecoyHis + 1.0 + snd + (fun (score, _) -> score) + (fun (score, _) -> score) + [| (-0.3, true) |] + |> Array.exactlyOne + expectFloatClose negativeBin -0.5 "the negative bin label uses the negative half-bandwidth branch" + + testCase "calculatePEPValues returns score-sorted decoy/total ratios" <| fun _ -> + let dataFreq = [| (2.0, 4.0, 1.0); (1.0, 2.0, 1.0) |] + let actual = + FDRControl.calculatePEPValues + (fun (_, total, _) -> total) + (fun (_, _, decoy) -> decoy) + (fun (score, _, _) -> score) + dataFreq + Expect.equal actual [ (1.0, 0.5); (2.0, 0.25) ] "PEP values are sorted ascending by score" + + testCase "logitTransformPepValues drops endpoint pep values and log10-logit-transforms the rest" <| fun _ -> + let scores, pepValues = + FDRControl.logitTransformPepValues + [| 1.0; 2.0; 3.0; 4.0 |] + [| 0.0; 0.5; 0.9; 1.0 |] + Expect.equal scores [| 2.0; 3.0 |] "endpoint PEP values are removed" + expectFloatClose pepValues.[0] 0.0 "the PEP 0.5 logit is zero" + expectFloatClose pepValues.[1] (log10 9.0) "the PEP 0.9 logit is log10(9)" + + testCase "initCalculateLin yields a monotonically usable PEP mapping on a separable target/decoy set" <| fun _ -> + let targets = [| for score in 1.0 .. 0.25 .. 8.0 -> (score, false) |] + let decoys = + Array.append + [| for score in -8.0 .. 0.25 .. -1.0 -> (score, true) |] + [| (1.1, true); (1.2, true); (2.1, true); (3.1, true) |] + let data = Array.append targets decoys + let msgs = ResizeArray() + let f = + FDRControl.initCalculateLin + msgs.Add + 0.5 + snd + (fun (score, _) -> score) + (fun (score, _) -> score) + data + let atOne = f 1.0 + let atEight = f 8.0 + let atFive = f 5.0 + Expect.isTrue (not (Double.IsNaN atOne) && not (Double.IsInfinity atOne)) "the PEP at score 1.0 is finite" + Expect.isTrue (not (Double.IsNaN atEight) && not (Double.IsInfinity atEight)) "the PEP at score 8.0 is finite" + Expect.isTrue (atEight <= atOne) "the PEP does not increase for the better score" + Expect.isTrue (0.0 <= atFive && atFive <= 1.0) "the PEP at score 5.0 is a probability" + Expect.equal msgs.Count 3 "the initializer emits the two setup traces and chosen bandwidth trace" + + testCase "getLogisticRegressionFunction fits a descending logistic mapping" <| fun _ -> + let f = + FDRControl.getLogisticRegressionFunction + (vector [| 1.0; 2.0; 3.0; 4.0; 5.0; 6.0 |]) + (vector [| 1.0; 1.0; 1.0; 0.0; 0.0; 0.0 |]) + 0.0001 + let atMiddle = f 3.5 + Expect.isTrue (0.0 <= atMiddle && atMiddle <= 1.0) "the logistic prediction is a probability" + Expect.isTrue (f 1.0 > f 6.0) "the logistic prediction is higher on the y=1 side" + ] + testList "LogisticRegression" [ testCase "calculateQValueLogReg yields a finite, broadly descending q-value function on well-separated data" <| fun _ -> let targets = diff --git a/tests/BioFSharp.Mz.Tests/ProteinInferenceTests.fs b/tests/BioFSharp.Mz.Tests/ProteinInferenceTests.fs index 14e880c..3a0730d 100644 --- a/tests/BioFSharp.Mz.Tests/ProteinInferenceTests.fs +++ b/tests/BioFSharp.Mz.Tests/ProteinInferenceTests.fs @@ -149,6 +149,158 @@ let tests = Expect.isSome (ProteinInference.isRNA rnaLine) "isRNA recognizes the mRNA feature" Expect.isNone (ProteinInference.isRNA geneLine) "isRNA rejects the gene feature" // The documented feature-string predicates classify GFF lines by feature type. + + testCase "createProteinModelInfoFromEntry reads the ID attribute as splice variant id" <| fun _ -> + let entry : GFFEntry = + { + Seqid = "chr1" + Source = "source" + Feature = "mRNA" + StartPos = 1 + EndPos = 10 + Score = 0.0 + Strand = '+' + Phase = 0 + Attributes = Map.ofList ["ID", ["Cre01.g000050.t1.1"]] + Supplement = [||] + } + let modelInfo = ProteinInference.createProteinModelInfoFromEntry 0 "locusA" entry + Expect.equal modelInfo.Id "Cre01.g000050.t1.1" "the ID attribute is used as the splice variant id" + Expect.equal modelInfo.Strand StrandDirection.Forward "a plus strand is mapped to Forward" + // GFF3's ID attribute identifies the splice variant that is matched against the protein sequence. + + testCase "createProteinModelInfoFromEntry tolerates unknown strand characters as Forward" <| fun _ -> + let mkEntry strand : GFFEntry = + { + Seqid = "chr1" + Source = "source" + Feature = "mRNA" + StartPos = 1 + EndPos = 10 + Score = 0.0 + Strand = strand + Phase = 0 + Attributes = Map.ofList ["ID", ["Cre01.g000050.t1.1"]] + Supplement = [||] + } + let unknownStrandInfo = + ProteinInference.createProteinModelInfoFromEntry 0 "locusA" (mkEntry '.') + let reverseStrandInfo = + ProteinInference.createProteinModelInfoFromEntry 0 "locusA" (mkEntry '-') + Expect.equal unknownStrandInfo.Strand StrandDirection.Forward "an unknown strand character falls back to Forward" + Expect.equal reverseStrandInfo.Strand StrandDirection.Reverse "a minus strand is mapped to Reverse" + // Unknown GFF strand characters use the documented Forward fallback while '-' remains Reverse. + + testCase "createProteinModelInfoFromEntry fails when the ID attribute is missing" <| fun _ -> + let entry : GFFEntry = + { + Seqid = "chr1" + Source = "source" + Feature = "mRNA" + StartPos = 1 + EndPos = 10 + Score = 0.0 + Strand = '+' + Phase = 0 + Attributes = Map.ofList ["Name", ["x"]] + Supplement = [||] + } + Expect.throws + (fun () -> ProteinInference.createProteinModelInfoFromEntry 0 "locusA" entry |> ignore) + "an entry without an ID attribute is rejected" + // The reader now requires ID and no longer accepts the legacy Name attribute. + + testCase "assignTranscriptsToGenes maps transcripts through the supplied parser" <| fun _ -> + let geneEntry : GFFEntry = + { + Seqid = "chr1" + Source = "source" + Feature = "gene" + StartPos = 1 + EndPos = 10 + Score = 0.0 + Strand = '+' + Phase = 0 + Attributes = Map.ofList ["ID", ["gene1"]] + Supplement = [||] + } + let rnaEntry : GFFEntry = + { + Seqid = "chr1" + Source = "source" + Feature = "mRNA" + StartPos = 1 + EndPos = 10 + Score = 0.0 + Strand = '+' + Phase = 0 + Attributes = Map.ofList ["ID", ["gene1.t1"]] + Supplement = [||] + } + let geneLine : GFFLine> = GFFEntryLine geneEntry + let rnaLine : GFFLine> = GFFEntryLine rnaEntry + let result = + ProteinInference.assignTranscriptsToGenes + (fun (s: string) -> Some (s + "!")) + [geneLine; rnaLine] + Expect.equal result.Count 1 "one transcript is assigned" + Expect.isTrue (result.ContainsKey "gene1.t1!") "the supplied parser determines the map key" + Expect.equal result.["gene1.t1!"].Id "gene1.t1" "the model info retains the original GFF3 ID" + // Transcript IDs are transformed by the supplied parser before entering the protein-model map. + + testCase "assignTranscriptsToGenes fails when the parser rejects an id" <| fun _ -> + let geneEntry : GFFEntry = + { + Seqid = "chr1" + Source = "source" + Feature = "gene" + StartPos = 1 + EndPos = 10 + Score = 0.0 + Strand = '+' + Phase = 0 + Attributes = Map.ofList ["ID", ["gene1"]] + Supplement = [||] + } + let rnaEntry : GFFEntry = + { + Seqid = "chr1" + Source = "source" + Feature = "mRNA" + StartPos = 1 + EndPos = 10 + Score = 0.0 + Strand = '+' + Phase = 0 + Attributes = Map.ofList ["ID", ["gene1.t1"]] + Supplement = [||] + } + let geneLine : GFFLine> = GFFEntryLine geneEntry + let rnaLine : GFFLine> = GFFEntryLine rnaEntry + Expect.throws + (fun () -> + ProteinInference.assignTranscriptsToGenes + (fun (_: string) -> None) + [geneLine; rnaLine] + |> ignore) + "a rejected transcript ID raises an exception" + // Every transcript must match the supplied parser or assignment fails with the documented error. + + testCase "PSMInput maps its score field to the ModelScore column" <| fun _ -> + let scoreField = + Microsoft.FSharp.Reflection.FSharpType.GetRecordFields typeof + |> Array.find (fun field -> field.Name = "Score") + let fieldAttributeType = typeof + let fieldAttribute = + scoreField.GetCustomAttributes(fieldAttributeType, false) + |> Array.exactlyOne + Expect.isNotNull fieldAttribute "the Score field has a FieldAttribute" + let columnIdentifier = + System.Reflection.CustomAttributeData.GetCustomAttributes(scoreField) + |> Seq.find (fun attribute -> attribute.AttributeType = fieldAttributeType) + |> fun attribute -> attribute.ConstructorArguments.[0].Value :?> string + Expect.equal columnIdentifier "ModelScore" "the Score field is read from the ModelScore column" + // FieldAttribute does not publicly expose its original string identifier, so the constructor metadata is inspected. ] testList "Inference" [ diff --git a/tests/BioFSharp.Mz.Tests/SearchDBTests.fs b/tests/BioFSharp.Mz.Tests/SearchDBTests.fs index 62466f2..d07cdf1 100644 --- a/tests/BioFSharp.Mz.Tests/SearchDBTests.fs +++ b/tests/BioFSharp.Mz.Tests/SearchDBTests.fs @@ -481,6 +481,45 @@ let tests = Expect.isTrue (List.contains "LLVR" plainSequences) "the configured minimum-length LLVR peptide is present" ) + testCase "getThreadSafePeptideLookUpFromFileBySequenceAndGMod retrieves the same entry as the mass-range lookup" <| fun _ -> + withTemporaryDirectory (fun directory -> + let fastaPath = Path.Combine(directory, "t8sg.fasta") + File.WriteAllText(fastaPath, ">sp|TESTPROT1|TEST\r\nMAGSTKLLVR\r\n") + + let sdbParams = + SearchDB.createSearchDbParams + "t8sg" + directory + fastaPath + id + (Digestion.Table.getProteaseBy "Trypsin") + 0 + 0 + 2000.0 + 4 + 20 + [] + SearchDB.MassMode.Monoisotopic + monoisotopicMass + [] + [SearchDB.Table.oxidation'Met'] + 1 + + use connection = SearchDB.connectOrCreateDB sdbParams + let byRange = SearchDB.getThreadSafePeptideLookUpFromFileBy connection sdbParams + let expected = + byRange 0.0 2000.0 + |> List.sortBy (fun result -> result.ModSequenceID) + |> List.head + let bySeq = SearchDB.getThreadSafePeptideLookUpFromFileBySequenceAndGMod connection sdbParams + let actual = bySeq expected.StringSequence expected.GlobalMod + + Expect.equal actual.PepSequenceID expected.PepSequenceID "the sequence-and-global-mod lookup returns the same peptide sequence" + Expect.equal actual.StringSequence expected.StringSequence "the sequence-and-global-mod lookup returns the same sequence" + Expect.equal actual.GlobalMod expected.GlobalMod "the sequence-and-global-mod lookup returns the same global modification" + expectWithin 0.000001 actual.Mass expected.Mass "the sequence-and-global-mod lookup returns the same mass" + ) + // PENDING: the function opens an already-open connection (getDBConnectionBy opens; the // function calls Open() again) and throws InvalidOperationException on every call. ptestCase "getProteinLookUpFromFileBy returns the protein for a peptide sequence ID" <| fun _ -> diff --git a/tests/BioFSharp.Mz.Tests/SparsePeakArrayTests.fs b/tests/BioFSharp.Mz.Tests/SparsePeakArrayTests.fs new file mode 100644 index 0000000..47174ce --- /dev/null +++ b/tests/BioFSharp.Mz.Tests/SparsePeakArrayTests.fs @@ -0,0 +1,51 @@ +module SparsePeakArrayTests + +open Expecto +open BioFSharp.Mz + +[] +let tests = + testList "SparsePeakArrayTests" [ + testCase "initMzToBinIdx and initBinIdxToMz are consistent for exact bin centers" <| fun _ -> + let width = 0.5 + let offset = 0.4 + let binIdx = SparsePeakArray.initMzToBinIdx width offset 100.3 + let mz = SparsePeakArray.initBinIdxToMz width offset binIdx + // Exact bin-center conversion must preserve the expected index and m/z. + Expect.equal binIdx 201 "100.3 maps to bin 201" + Expect.floatClose Accuracy.high mz 100.3 "bin 201 maps back to m/z 100.3" + + testCase "peaksToNearestBinVector sums intensities of peaks falling in the same bin" <| fun _ -> + let peaks = [| Peak(10.2, 5.0); Peak(10.4, 3.0); Peak(12.7, 2.0) |] + let result = SparsePeakArray.peaksToNearestBinVector 1.0 0.0 0.0 1000.0 peaks + let bin10 = SparsePeakArray.initMzToBinIdx 1.0 0.0 10.2 + let bin12 = SparsePeakArray.initMzToBinIdx 1.0 0.0 12.7 + // Peaks sharing a bin are aggregated by summing their intensities. + Expect.equal result.Data.Count 2 "only two occupied bins are present" + Expect.isTrue (result.Data.ContainsKey bin10) "bin 10 is present" + Expect.isTrue (result.Data.ContainsKey bin12) "bin 12 is present" + Expect.floatClose Accuracy.high result.Data.[bin10] 8.0 "bin 10 contains the summed intensity" + Expect.floatClose Accuracy.high result.Data.[bin12] 2.0 "bin 12 contains its peak intensity" + + testCase "peaksToNearestBinVector excludes peaks at or beyond the mass borders" <| fun _ -> + let peaks = [| Peak(10.0, 1.0); Peak(20.0, 2.0); Peak(15.0, 4.0) |] + let result = SparsePeakArray.peaksToNearestBinVector 1.0 0.0 10.0 20.0 peaks + let bin15 = SparsePeakArray.initMzToBinIdx 1.0 0.0 15.0 + // The lower and upper mass borders are strict exclusion boundaries. + Expect.equal result.Data.Count 1 "only the interior peak's bin is present" + Expect.isTrue (result.Data.ContainsKey bin15) "the interior peak's bin is present" + Expect.floatClose Accuracy.high result.Data.[bin15] 4.0 "the interior peak is retained" + + testCase "dot multiplies matching bins and ignores disjoint ones" <| fun _ -> + let x = + SparsePeakArray.peaksToNearestBinVector 1.0 0.0 0.0 1000.0 [| Peak(10.2, 2.0); Peak(12.7, 3.0) |] + let y = + SparsePeakArray.peaksToNearestBinVector 1.0 0.0 0.0 1000.0 [| Peak(10.3, 4.0); Peak(30.0, 7.0) |] + let disjointX = + SparsePeakArray.peaksToNearestBinVector 1.0 0.0 0.0 1000.0 [| Peak(100.0, 2.0) |] + let disjointY = + SparsePeakArray.peaksToNearestBinVector 1.0 0.0 0.0 1000.0 [| Peak(200.0, 4.0) |] + // Only bin 10 is shared, so the dot product is 2.0 * 4.0. + Expect.floatClose Accuracy.high (SparsePeakArray.dot x y) 8.0 "matching bins are multiplied and summed" + Expect.floatClose Accuracy.high (SparsePeakArray.dot disjointX disjointY) 0.0 "disjoint bins contribute nothing" + ]