BiobaseBlast 0.3.1.0 → 0.3.3.0

raw patch · 10 files changed

+135/−63 lines, 10 filesdep ~BiobaseTypesdep ~BiobaseXNAdep ~PrimitiveArraynew-uploaderPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: BiobaseTypes, BiobaseXNA, PrimitiveArray, SciBaseTypes

API changes (from Hackage documentation)

Files

• Biobase/BLAST/Import.hs +3/−3
• Biobase/BLAST/Types.hs +28/−30
• Biobase/SubstMatrix.hs +18/−9
• Biobase/SubstMatrix/Embedded.hs +3/−2
• Biobase/SubstMatrix/Import.hs +21/−8
• Biobase/SubstMatrix/Statistics.hs +41/−0
• Biobase/SubstMatrix/Types.hs +6/−3
• BiobaseBlast.cabal +7/−6
• README.md +1/−2
• changelog.md +7/−0

@@ -36,7 +36,7 @@        bs <- B.readFile filePath        let json = parseJSONBlastCmd bs        return json-     else fail "# JSON2 blast file \"%s\" does not exist\n" filePath+     else error "# JSON2 blast file \"%s\" does not exist\n" filePath  parseJSONBlastCmd :: B.ByteString -> Either String BlastCmdJSON2 parseJSONBlastCmd bs = A.eitherDecode bs :: Either String BlastCmdJSON2@@ -51,7 +51,7 @@   blastFileExists <- doesFileExist filePath   if blastFileExists      then parseTabularBlasts <$> B.readFile filePath-     else fail "# tabular blast file \"%s\" does not exist\n" filePath+     else error "# tabular blast file \"%s\" does not exist\n" filePath  -- | reads and parses tabular HTTP Blast result from provided filePath blastHTTPFromFile :: String -> IO [BlastTabularResult]@@ -60,7 +60,7 @@   blastFileExists <- doesFileExist filePath   if blastFileExists       then parseTabularHTTPBlasts <$> B.readFile filePath-      else fail "# tabular blast file \"%s\" does not exist\n" filePath+      else error "# tabular blast file \"%s\" does not exist\n" filePath  -- | Read a lazy bytestring and stream out a lsit of @BlastTabularResult@'s. -- In case, there is a parse error "late" in the file, we might have

@@ -26,7 +26,8 @@ import Data.Aeson.Types import qualified Data.HashMap.Strict as HM import qualified Data.Sequence as DS-import Control.Lens+import Control.Lens(makeLenses)+import Data.HashMap.Strict (elems)  -- | Turn all keys in a JSON object to lowercase. jsonLower :: Value -> Value@@ -204,15 +205,37 @@ data BlastProgram = BlastX | BlastP | BlastN   deriving (Show, Eq) +makeLenses ''BlastTabularResult+makeLenses ''BlastTabularHit+makeLenses ''SearchTarget+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''SearchTarget+makeLenses ''BlastProgram+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastProgram+makeLenses ''SearchStat+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''SearchStat+makeLenses ''HitDescription+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''HitDescription+makeLenses ''Hit+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''Hit+makeLenses ''Hsp+makeLenses ''Search+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''Search+makeLenses ''BlastJSONResult+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastJSONResult+makeLenses ''Params+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''Params+makeLenses ''BlastReport+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastReport+makeLenses ''BlastOutput2+deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastOutput2 makeLenses ''BlastCmdJSON2 --deriveJSON defaultOptions{fieldLabelModifier = (map toLower) . drop 1} ''BlastJSON2 instance FromJSON BlastCmdJSON2 where-  parseJSON (Object v) =-    BlastCmdJSON2 <$> (v .: "BlastOutput2")+  parseJSON (Object v) = BlastCmdJSON2 <$> (v .: "BlastOutput2")   parseJSON _ = mzero instance ToJSON BlastCmdJSON2 where   toJSON (BlastCmdJSON2 _blastoutput2) =-        object [ "BlastOutput2"  Data.Aeson.Types..= _blastoutput2 ]+        object [ "BlastOutput2" .= _blastoutput2 ]  makeLenses ''BlastJSON2 --deriveJSON defaultOptions{fieldLabelModifier = (map toLower) . drop 1} ''BlastJSON2@@ -221,29 +244,4 @@   parseJSON _ = mzero  instance ToJSON BlastJSON2 where-  toJSON (BlastJSON2 _blastoutput2) = object [ "BlastOutput2"  Data.Aeson.Types..= _blastoutput2 ]---deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastJSON2-makeLenses ''BlastOutput2-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastOutput2-makeLenses ''BlastReport-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastReport-makeLenses ''Params-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''Params-makeLenses ''BlastJSONResult-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastJSONResult-makeLenses ''Search-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''Search-makeLenses ''Hit-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''Hit-makeLenses ''Hsp---deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''Hsp-makeLenses ''HitDescription-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''HitDescription-makeLenses ''SearchStat-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''SearchStat-makeLenses ''BlastTabularResult-makeLenses ''BlastTabularHit-makeLenses ''SearchTarget-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''SearchTarget-makeLenses ''BlastProgram-deriveJSON defaultOptions{fieldLabelModifier = drop 1} ''BlastProgram+  toJSON (BlastJSON2 _blastoutput2) = object [ "BlastOutput2"  .= _blastoutput2 ]

@@ -10,7 +10,9 @@ import           Data.List (maximumBy,find) import           Data.Serialize (Serialize) import           Data.Vector.Unboxed.Deriving+import           Debug.Trace (traceShow) import           GHC.Generics (Generic)+import           GHC.Real(Ratio(..)) import           Numeric.Log import qualified Data.Map.Strict as M import qualified Data.Vector.Unboxed as VU@@ -36,6 +38,7 @@  import           Biobase.SubstMatrix.Embedded import           Biobase.SubstMatrix.Import+import           Biobase.SubstMatrix.Statistics import           Biobase.SubstMatrix.Types  @@ -46,28 +49,32 @@ -- The resulting @AA@ are tagged with the name type from the @DNA n@. -- -- TODO Definitely use the correct upper bound constants here!+--+-- TODO the amino-acid range produces wrong reads from @m@. We should have this in ExceptT to make+-- sure we do not read from undefined memory.  mkANuc3SubstMat   ∷ TranslationTable (Letter DNA n) (Letter AA n)-  → AASubstMat t (DiscLogOdds k) n+  → AASubstMat t (DiscLogOdds k) n n   → ANuc3SubstMat t (Letter AA n, (DiscLogOdds k)) n n-mkANuc3SubstMat tbl (AASubstMat m)+mkANuc3SubstMat tbl (AASubstMat m l)   = ANuc3SubstMat   $ fromAssocs (ZZ:..LtLetter AA.Undef:..LtLetter DNA.N:..LtLetter DNA.N:..LtLetter DNA.N) (AA.Undef, DiscLogOdds . Discretized $ -999)     [ ( (Z:.a:.u:.v:.w)-      , (t, m!(Z:.a:.t))+      , (t, x)       )     | a <- VU.toList aaRange     , u <- [DNA.A .. DNA.N], v <- [DNA.A .. DNA.N], w <- [DNA.A .. DNA.N]     , let b = Codon u v w     , let t = translate tbl b+    , x <- maybe [] (:[]) $ m!?(Z:.a:.t)     ]  -- | This function does the following: -- 1. check if @fname@ is a file, and if so try to load it. -- 2. if not, check if @fname@ happens to be the name of one of the known @PAM/BLOSUM@ tables. -fromFileOrCached ∷ (MonadIO m, MonadError String m) ⇒ FilePath → m (AASubstMat t (DiscLogOdds Unknown) a)+fromFileOrCached ∷ (MonadIO m, MonadError String m) ⇒ FilePath → m (AASubstMat t (DiscLogOdds (1 :% 1)) a b) fromFileOrCached fname = do   dfe ← liftIO $ doesFileExist fname   if | fname == "list" → do@@ -79,13 +86,15 @@  -- | Turn log-odds into log-probabilities. Normalizes over the whole set of -- values in the matrix.+--+-- TODO check if lambda should indeed be 1.0  mkProbabilityMatrix-  ∷ Double-  → AASubstMat t (DiscLogOdds k) n-  → AASubstMat t (Log (Probability NotNormalized Double)) n-mkProbabilityMatrix invScale (AASubstMat dlo) = AASubstMat $ PA.map (/nrm) $ dbl-  where dbl = PA.map (\(DiscLogOdds (Discretized k)) → stateLogProbability (negate invScale) $ fromIntegral @Int @Double k) dlo+  :: Double+  -> AASubstMat t (DiscLogOdds k) n m+  -> AASubstMat t (Log (Probability NotNormalized Double)) n m+mkProbabilityMatrix invScale (AASubstMat dlo l) = AASubstMat (mapArray (/nrm) $ dbl) 1.0+  where dbl = mapArray (\(DiscLogOdds (Discretized k)) → stateLogProbability (negate invScale) $ fromIntegral @Int @Double k) dlo         nrm = maximum . Prelude.map snd $ PA.assocs dbl  

@@ -5,6 +5,7 @@ import Data.ByteString (ByteString) import Data.FileEmbed import Control.Arrow (second)+import GHC.Real(Ratio(..))  import Numeric.Discretized import Statistics.Odds@@ -15,9 +16,9 @@   embeddedPamBlosumFiles ∷ [(FilePath,ByteString)]-embeddedPamBlosumFiles = $(embedDir "sources/PamBlosum")+embeddedPamBlosumFiles = $(makeRelativeToProject "sources/PamBlosum" >>= embedDir) -embeddedPamBlosum ∷ [(FilePath,AASubstMat t (DiscLogOdds Unknown) a)]+embeddedPamBlosum ∷ [(FilePath,AASubstMat t (DiscLogOdds (1:%1)) a b)] embeddedPamBlosum = either error id . runExcept                   . mapM (\(k,v) → fromByteString v >>= \mv → return (k,mv))                   $ embeddedPamBlosumFiles

@@ -20,6 +20,7 @@ import           Statistics.Odds  import           Biobase.SubstMatrix.Types+import           Biobase.SubstMatrix.Statistics   @@ -27,19 +28,31 @@ -- -- TODO the parser is fragile, since it uses @read@. This should be fixed. -fromByteString ∷ (MonadError String m) ⇒ ByteString → m (AASubstMat t (DiscLogOdds k) a)-fromByteString bs = do-  let (x:xs) = dropWhile (("#"==).take 1) . lines $ unpack bs-  let cs = map head . words $ x -- should give us the characters encoding an amino acid-  let ss = map (map (DiscLogOdds . Discretized) . map read . drop 1 . words) $ xs-  let xs = [ ((Z:.charAA k1:.charAA k2),z)+fromByteString+  :: ( MonadError String m+     , Real (Discretized k)+     )+  => ByteString+  -> m (AASubstMat t (DiscLogOdds k) a b)+fromByteString bstring = do+  let (b:bs) = dropWhile (("#"==).take 1) . lines $ unpack bstring+      cs = map head . words $ b -- should give us the characters encoding an amino acid+      ss = map (map (DiscLogOdds . Discretized) . map read . drop 1 . words) $ bs+      xs = [ ((Z:.charAA k1:.charAA k2),z)            | (k1,s) <- zip cs ss            , (k2,z) <- zip cs s            ]-  return . AASubstMat $ fromAssocs (ZZ:..LtLetter AA.Z:..LtLetter AA.Z) (DiscLogOdds . Discretized $ -999) xs+      as = fromAssocs (ZZ:..LtLetter AA.Z:..LtLetter AA.Z) (DiscLogOdds . Discretized $ -999) xs+      mat = AASubstMat as (estimateLambda mat)+  return mat  -- | Import substitution matrix from file. -fromFile ∷ (MonadIO m, MonadError String m) ⇒ FilePath → m (AASubstMat t (DiscLogOdds k) a)+fromFile+  :: ( MonadIO m, MonadError String m+     , Real (Discretized k)+     )+  => FilePath+  -> m (AASubstMat t (DiscLogOdds k) a b) fromFile fname = liftIO (BS.readFile fname) >>= fromByteString 

@@ -0,0 +1,41 @@++-- | This module provides statistics needed for substitution matrices. It is a+-- very modest attempt to replicate some of the Blast statistics.++module Biobase.SubstMatrix.Statistics where++import Data.Vector.Unboxed (Unbox)+import Data.Vector.Fusion.Util+import Data.Vector.Fusion.Stream.Monadic as SM+import Debug.Trace++import Data.PrimitiveArray as PA+import           Biobase.Primary.Letter++import Biobase.SubstMatrix.Types++++-- | estimate Blast lambda.+--+-- TODO use ExceptT++estimateLambda+  :: (Unbox s, Num s, Real s)+  => AASubstMat t s a b -> Double+{-# Inlinable estimateLambda #-}+estimateLambda (AASubstMat mat _) = go 1000 1 2 0 where+  go count lambda high low+    | count <= 0 = error "no convergence?!"+    | (high-low) <= 0.001 = lambda+    | s >  1              = go (count-1) ((lambda+low)/2)  lambda low+    | s <= 1              = go (count-1) ((lambda+high)/2) high   lambda+    where+      -- get the rows and columns, needed to get the probs for each row/column right.+      (ZZ:..r':..c') = PA.upperBound mat+      r = fromRational $ toRational $ size r'+      c = fromRational $ toRational $ size c'+      -- sum of all scores+      s = unId . SM.foldl' (+) 0 $ SM.map eachElem $ PA.assocsS mat+      eachElem (Z:.i:.j, z) = (1/r) * (1/c) * exp (lambda * (fromRational $ toRational z))+

@@ -27,7 +27,10 @@  -- An amino-acid substitution matrix. Tagged with the type of scoring used. -newtype AASubstMat t s a = AASubstMat { _aaSubstMat :: Unboxed (Z:.Letter AA a:.Letter AA a) s }+data AASubstMat t s a b = AASubstMat+  { _aaSubstMat     :: !(Unboxed (Z:.Letter AA a:.Letter AA b) s)+  , _aaSubstLambda  :: Double+  }   deriving (Generic,Eq,Read,Show) makeLenses ''AASubstMat @@ -40,11 +43,11 @@  -- | @PAM@ matrices are similarity matrices. -type SubstPAM = AASubstMat Similarity (DiscLogOdds Unknown)+type SubstPAM = AASubstMat Similarity (DiscLogOdds (RTyId 1))  -- | @BLOSUM@ matrices are distance matrices. -type SubstBLOSUM = AASubstMat Distance (DiscLogOdds Unknown)+type SubstBLOSUM = AASubstMat Distance (DiscLogOdds (RTyId 1))  -- | Substitution matrix from amino acids to nucleotide triplets. 

@@ -1,5 +1,5 @@ name:           BiobaseBlast-version:        0.3.1.0+version:        0.3.3.0 author:         Christian Hoener zu Siederdissen, Florian Eggenhofer maintainer:     choener@bioinf.uni-leipzig.de homepage:       https://github.com/choener/BiobaseBlast@@ -11,7 +11,7 @@ build-type:     Simple stability:      experimental cabal-version:  >= 1.10.0-tested-with:    GHC == 8.4.4+tested-with:    GHC == 8.8, GHC == 8.10, GHC == 9.0 synopsis:       BLAST-related tools description:                 This library contains BLAST-related functionality:@@ -151,10 +151,10 @@                , vector-th-unbox  >= 0.2                --                , BiobaseENA       == 0.0.0.*-               , BiobaseTypes     == 0.2.0.*-               , BiobaseXNA       == 0.11.0.*-               , PrimitiveArray   == 0.9.1.*-               , SciBaseTypes     == 0.1.0.*+               , BiobaseTypes     == 0.2.1.*+               , BiobaseXNA       == 0.11.1.*+               , PrimitiveArray   == 0.10.1.*+               , SciBaseTypes     == 0.1.1.*    default-language:     Haskell2010@@ -178,6 +178,7 @@     Biobase.SubstMatrix.Embedded     Biobase.SubstMatrix.Hints     Biobase.SubstMatrix.Import+    Biobase.SubstMatrix.Statistics     Biobase.SubstMatrix.Types    ghc-options:

@@ -1,6 +1,5 @@-[![Build Status](https://travis-ci.org/choener/BiobaseBlast.svg?branch=master)](https://travis-ci.org/choener/BiobaseBlast)+# BiobaseBlast [![Hackage](https://img.shields.io/hackage/v/BiobaseBlast.svg)](https://hackage.haskell.org/package/BiobaseBlast) ![github action: CI](https://github.com/choener/BiobaseBlast/actions/workflows/ci.yml/badge.svg) -# BiobaseBlast  This library contains BLAST-related functionality. 

@@ -1,3 +1,10 @@+0.3.3.0+-------++- Compatibility with GHC 9+- Switched testing to github actions++ 0.3.1.0 -------