packages feed

BiobaseBlast 0.0.0.1 → 0.2.0.0

raw patch · 10 files changed

+443/−66 lines, 10 filesdep +BiobaseBlastdep +BiobaseTypesdep +PrimitiveArraydep −arraydep ~BiobaseXNAdep ~base

Dependencies added: BiobaseBlast, BiobaseTypes, PrimitiveArray, aeson, attoparsec, binary, bytestring, cereal, deepseq, directory, filepath, split, tasty, tasty-quickcheck, tasty-silver, tasty-th, vector, vector-th-unbox

Dependencies removed: array

Dependency ranges changed: BiobaseXNA, base

Files

+ Biobase/BLAST.hs view
@@ -0,0 +1,11 @@+-- | Types and functions for NCBI BLAST++--++module Biobase.BLAST+  ( module Biobase.BLAST.Types+  , module Biobase.BLAST.Import+  ) where++import Biobase.BLAST.Import (blastFromFile)+import Biobase.BLAST.Types+
+ Biobase/BLAST/Import.hs view
@@ -0,0 +1,130 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE LambdaCase #-}++-- | Parses NCBI BLAST+ tabular output++module Biobase.BLAST.Import where++import Prelude hiding (takeWhile)+import Data.Attoparsec.ByteString.Char8 hiding (isSpace)+import qualified Data.Attoparsec.ByteString.Lazy as L+import qualified Data.ByteString.Char8 as C+import qualified Data.ByteString.Builder as S+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.Vector as V+import System.Directory+import Data.Char+import Control.Monad+import Debug.Trace+import Text.Printf+import Biobase.BLAST.Types++-- | reads and parses tabular Blast result from provided filePath+blastFromFile :: String -> IO [BlastTabularResult]+blastFromFile filePath = do+  printf "# reading tabular blast input from file %s\n" filePath+  blastFileExists <- doesFileExist filePath+  if blastFileExists+     then parseTabularBlasts <$> B.readFile filePath+     else fail "# 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+-- already streamed out some (or many!) of these results.++parseTabularBlasts :: B.ByteString -> [BlastTabularResult]+parseTabularBlasts = go+  where go xs = case L.parse genParseTabularBlast xs of+          L.Fail remainingInput ctxts err  -> error $ "parseTabularBlasts failed! " ++ err ++ " ctxt: " ++ show ctxts ++ " head of remaining input: " ++ B.unpack (B.take 1000 remainingInput)+          L.Done remainingInput btr+            | B.null remainingInput  -> [btr]+            | otherwise              -> btr : go remainingInput++genParseBlastProgram :: Parser BlastProgram+genParseBlastProgram = do+  choice [string "# BLAST",string "# blast"]+  (toLower <$> anyChar) >>= return . \case +    'x' -> BlastX+    'p' -> BlastP+    'n' -> BlastN++genParseTabularBlast :: Parser BlastTabularResult+genParseTabularBlast = do+  --choice [string "# BLAST",string "# blast"]+  --_blastProgram <- choice [string "p",string "P",string "x",string "X"] <?> "Program"+  _blastProgram <- genParseBlastProgram <?> "Program"+  many1 (notChar '\n')+  endOfLine+  string "# Query: " <?> "Query"+  --_blastQueryId <- many1 (notChar ' ') <* skipWhile (not (string "\n")) <?> "QueryId"+  _blastQueryId <- takeWhile (not . isSpace) <* manyTill anyChar endOfLine <?> "QueryId"+  --_blastQueryName <- many' (notChar '\n')  <?> "QueryName"+  string "# Database: " <?> "Database"+  _blastDatabase <- many1 (notChar '\n') <?> "Db"+  string "\n# " <?> "header linebreak"+  --fields line+  skipMany (try genParseFieldLine) <?> "Fields"+  _blastHitNumber <- decimal  <?> "Hit number"+  string " hits found\n" <?> "hits found"+  _tabularHit <- count  _blastHitNumber (try genParseBlastTabularHit)  <?> "Tabular hit"+  return $ BlastTabularResult _blastProgram (toLB _blastQueryId) (B.pack _blastDatabase) _blastHitNumber (V.fromList _tabularHit)++genParseFieldLine :: Parser ()+genParseFieldLine = do+  string "Fields:"+  skipMany (notChar '\n')+  string "\n# "+  return ()++genParseBlastTabularHit :: Parser BlastTabularHit+genParseBlastTabularHit = do+  --start <- peekChar'+  --when (start == '#') (fail "irgendwas")+  --_queryId <- takeWhile (\c ->  c /= '#' || c /= '\t')  <?> "hit qid"+  _queryId <- takeWhile1 ((/=9) . ord) <?> "hit qid"+  char '\t'+  _subjectId <- takeWhile1 ((/=9) . ord) <?> "hit sid"+  char '\t'+  _seqIdentity <- double <?> "hit seqid"+  char '\t'+  _alignmentLength <- decimal  <?> "hit sid"+  char '\t'+  _misMatches <- decimal <?> "hit mmatch"+  char '\t'+  _gapOpenScore <- decimal <?> "hit gopen"+  char '\t'+  _queryStart <- decimal <?> "hit qstart"+  char '\t'+  _queryEnd <- decimal  <?> "hit qend"+  char '\t'+  _hitSeqStart <- decimal  <?> "hit sstart"+  char '\t'+  _hitSeqEnd <- decimal <?> "hit send"+  char '\t'+  _eValue <- double <?> "hit eval"+  char '\t'+  _bitScore <- double <?> "hit bs"+  char '\t'+  _subjectFrame <- decimal <?> "hit sF"+  char '\t'+  _querySeq <- takeWhile1 ((/=9) . ord) <?> "hit qseq" -- 9 == '\t'+  char '\t'+  _subjectSeq <- takeWhile1 ((/=10) . ord) <?> "hit subSeq" -- 10 == '\n'+  char '\n'+  return $ BlastTabularHit (B.fromStrict _queryId) (B.fromStrict _subjectId) _seqIdentity _alignmentLength _misMatches _gapOpenScore _queryStart _queryEnd _hitSeqStart _hitSeqEnd _eValue _bitScore _subjectFrame (B.fromStrict _querySeq) (B.fromStrict _subjectSeq)+  +--IUPAC amino acid with gap+--aminoacidLetters :: Char -> Bool+aminoacidLetters = inClass "ARNDCQEGHILMFPSTWYVBZX-"++--IUPAC nucleic acid characters with gap+--nucleotideLetters :: Char -> Bool+nucleotideLetters = inClass "AGTCURYSWKMBDHVN-."++--IUPAC nucleic acid characters with gap+--bioLetters :: Char -> Bool+bioLetters = inClass "ABCDEFGHIJKLMNOPQRSTUVWXYZ.-"+++toLB :: C.ByteString -> B.ByteString+toLB = S.toLazyByteString . S.byteString
+ Biobase/BLAST/Types.hs view
@@ -0,0 +1,49 @@++-- | Encoding of tabular NCBI BLAST+ output++module Biobase.BLAST.Types where++import Prelude hiding (takeWhile)+import Data.Attoparsec.ByteString.Char8 hiding (isSpace)+import qualified Data.Attoparsec.ByteString.Lazy as L+import qualified Data.ByteString.Char8 as C+import qualified Data.ByteString.Builder as S+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.Vector as V+import System.Directory+import Data.Char+import Control.Monad+import Debug.Trace+import Text.Printf++data BlastTabularResult = BlastTabularResult+  { blastProgram :: !BlastProgram,+    blastQueryId :: !B.ByteString,+--    blastQueryName :: !B.ByteString,+    blastDatabase :: !B.ByteString,+    blastHitNumber :: !Int,+    hitLines :: !(V.Vector BlastTabularHit)+  }+  deriving (Show, Eq)++data BlastProgram = BlastX | BlastP | BlastN+  deriving (Show, Eq)++data BlastTabularHit = BlastTabularHit+  { queryId :: !B.ByteString,+    subjectId ::  !B.ByteString,+    seqIdentity :: !Double,+    alignmentLength :: !Int,+    misMatches :: !Int,+    gapOpenScore :: !Int,+    queryStart :: !Int,+    queryEnd :: !Int,+    hitSeqStart :: !Int,+    hitSeqEnd :: !Int,+    eValue :: !Double,+    bitScore :: !Double,+    subjectFrame :: !Int,+    querySeq  :: !B.ByteString,+    subjectSeq  :: !B.ByteString+  }+  deriving (Show, Eq)
Biobase/SubstMatrix.hs view
@@ -1,58 +1,131 @@  module Biobase.SubstMatrix where -import qualified Data.Array.IArray  as A-import qualified Data.Array.Unboxed as U-import qualified Data.Map as M-import Data.Char (toLower)+import           Control.DeepSeq (NFData(..))+import           Data.Aeson (FromJSON,ToJSON)+import           Data.Binary (Binary)+import           Data.Serialize (Serialize)+import           Data.Vector.Unboxed.Deriving+import           GHC.Generics (Generic)+import qualified Data.Map.Strict as M+import qualified Data.Vector.Unboxed as VU -import Biobase.Primary-import Biobase.Codon+import           Biobase.Primary.AA (AA,aaRange)+import           Biobase.Primary.Letter+import           Biobase.Primary.Nuc.DNA (DNA)+import           Biobase.Primary.Trans (dnaAAmap)+import           Biobase.Types.Odds+import           Data.PrimitiveArray+import qualified Biobase.Primary.AA as AA+import qualified Biobase.Primary.Nuc.DNA as D   --- | Substitution table for one amino acid with another+-- | Denotes that we are dealing with a similarity score. Higher is more+-- similar. -type SubstMatrix = U.Array (Char,Char) Int+data Similarity --- | Substitution from three DNA nucleotides to an amino acid.+-- | Denotes that we are dealing with a distance score. Lower is more+-- similar. -type Nuc3SubstMatrix = U.Array (Nuc,Nuc,Nuc,Char) Int+data Distance -type Nuc2SubstMatrix = U.Array (Nuc,Nuc,Char) Int+-- An amino-acid substitution matrix. Tagged with the type of scoring used. -type Nuc1SubstMatrix = U.Array (Nuc,Char) Int+newtype AASubstMat t = AASubstMat { aaSubstMat :: Unboxed (Z:.Letter AA:.Letter AA) DLO }+  deriving (Generic,Eq,Read,Show) +instance Binary    (AASubstMat t)+instance Serialize (AASubstMat t)+--instance FromJSON  (AASubstMat t)+--instance ToJSON    (AASubstMat t) +instance NFData (AASubstMat t) --- *+-- | @PAM@ matrices are similarity matrices. -mkNuc3SubstMatrix :: SubstMatrix -> Nuc3SubstMatrix-mkNuc3SubstMatrix mat = A.accumArray (\_ z -> z) (-10) ((nN,nN,nN,'*'),(nT,nT,nT,'Z'))-  [ ( (a,b,c,k), case l of Just l' -> mat A.! (l',k) ; Nothing -> -10 )-  | a<-acgt, b<-acgt, c<-acgt-  , k<-['*' .. 'Z']-  , let abc = map (toLower . fromNuc) [a,b,c]-  , let l = M.lookup abc codonTable-  ]+type SubstPAM = AASubstMat Similarity -mkNuc2SubstMatrix :: (Int -> Int -> Int) -> (Int -> Int) -> SubstMatrix -> Nuc2SubstMatrix-mkNuc2SubstMatrix f g mat = A.accumArray f (-10) ((nN,nN,'*'),(nT,nT,'Z'))-  [ ( (x,y,k), case l of Just l' -> mat A.! (l',k) ; Nothing -> -10 )-  | a<-acgt, b<-acgt, c<-acgt-  , k<-['*' .. 'Z']-  , (x,y) <- [ (a,b), (a,c), (b,c) ]-  , let abc = map (toLower . fromNuc) [a,b,c]-  , let l = M.lookup abc codonTable-  ]+-- | @BLOSUM@ matrices are distance matrices. -mkNuc1SubstMatrix :: (Int -> Int -> Int) -> (Int -> Int) -> SubstMatrix -> Nuc1SubstMatrix-mkNuc1SubstMatrix f g mat = A.accumArray f (-10) ((nN,'*'),(nT,'Z'))-  [ ( (x,k), case l of Just l' -> mat A.! (l',k) ; Nothing -> -10 )-  | a<-acgt, b<-acgt, c<-acgt-  , k<-['*' .. 'Z']-  , x <- [a,b,c]-  , let abc = map (toLower . fromNuc) [a,b,c]-  , let l = M.lookup abc codonTable+type SubstBLOSUM = AASubstMat Distance++-- | Substitution matrix from amino acids to nucleotide triplets.++newtype ANuc3SubstMat t = ANuc3SubstMat { anuc3SubstMat :: Unboxed (Z:.Letter AA:.Letter DNA:.Letter DNA:.Letter DNA) DLO }+  deriving (Generic,Eq,Read,Show)++instance Binary    (ANuc3SubstMat t)+instance Serialize (ANuc3SubstMat t)+--instance FromJSON  (ANuc3SubstMat t)+--instance ToJSON    (ANuc3SubstMat t)++instance NFData (ANuc3SubstMat t)++-- | Substitution matrix from amino acids to degenerate nucleotide+-- 2-tuples. The third nucleotide letter is missing.++newtype ANuc2SubstMat t = ANuc2SubstMat { anuc2SubstMat :: Unboxed (Z:.Letter AA:.Letter DNA:.Letter DNA) DLO }+  deriving (Generic,Eq,Read,Show)++instance Binary    (ANuc2SubstMat t)+instance Serialize (ANuc2SubstMat t)+--instance FromJSON  (ANuc2SubstMat t)+--instance ToJSON    (ANuc2SubstMat t)++instance NFData (ANuc2SubstMat t)++-- | Substitution matrix from amino acids to degenerate nucleotide+-- 1-tuples. Two out of three nucleotides in a triplet are missing.++newtype ANuc1SubstMat t = ANuc1SubstMat { anuc1SubstMat :: Unboxed (Z:.Letter AA:.Letter DNA) DLO }+  deriving (Generic,Eq,Read,Show)++instance Binary    (ANuc1SubstMat t)+instance Serialize (ANuc1SubstMat t)+--instance FromJSON  (ANuc1SubstMat t)+--instance ToJSON    (ANuc1SubstMat t)++instance NFData (ANuc1SubstMat t)++-- | The usual substitution matrix, but here with a codon and an amino acid+-- to be compared.++mkANuc3SubstMat :: AASubstMat t -> ANuc3SubstMat t+mkANuc3SubstMat (AASubstMat m) = ANuc3SubstMat $ fromAssocs (Z:. AA.Stop :. D.A:.D.A:.D.A) (Z:. AA.Z :. D.N:.D.N:.D.N) (DLO $ -999)+  [ ( (Z:.a:.u:.v:.w) , maybe (DLO $ -999) (\b -> m!(Z:.a:.b)) $ M.lookup uvw dnaAAmap)+  | a <- aaRange+  , u <- [D.A .. D.N], v <- [D.A .. D.N], w <- [D.A .. D.N]+  , let uvw = VU.fromList [u,v,w]   ]++-- | Create a 2-tuple to amino acid substitution matrix. Here, @f@ combines+-- all to entries that have the same 2-tuple index.++mkANuc2SubstMat :: (DLO -> DLO -> DLO) -> AASubstMat t -> ANuc2SubstMat t+mkANuc2SubstMat f (AASubstMat m) = ANuc2SubstMat $ fromAssocs (Z:. AA.Stop :. D.A:.D.A) (Z:. AA.Z :. D.N:.D.N) (DLO $ -999)+  . M.assocs+  . M.fromListWith f+  $ [ ((Z:.a:.x:.y), maybe (DLO $ -999) (\k -> m!(Z:.a:.k)) $ M.lookup uvw dnaAAmap)+    | a <- aaRange+    , u <- [D.A .. D.N], v <- [D.A .. D.N], w <- [D.A .. D.N]+    , (x,y) <- [ (u,v), (u,w), (v,w) ]+    , let uvw = VU.fromList [u,v,w]+    ]++-- | The most degenerate case, where just a single nucleotide remains in+-- the amino-acid / nucleotide substitution. Again, @f@ combines different+-- entries.++mkANuc1SubstMat :: (DLO -> DLO -> DLO) -> AASubstMat t -> ANuc1SubstMat t+mkANuc1SubstMat f (AASubstMat m) = ANuc1SubstMat $ fromAssocs (Z:. AA.Stop :. D.A) (Z:. AA.Z :. D.N) (DLO $ -999)+  . M.assocs+  . M.fromListWith f+  $ [ ((Z:.a:.x), maybe (DLO $ -999) (\k -> m!(Z:.a:.k)) $ M.lookup uvw dnaAAmap)+    | a <- aaRange+    , u <- [D.A .. D.N], v <- [D.A .. D.N], w <- [D.A .. D.N]+    , x <- [u,v,w]+    , let uvw = VU.fromList [u,v,w]+    ]
Biobase/SubstMatrix/Import.hs view
@@ -1,17 +1,28 @@  module Biobase.SubstMatrix.Import where -import qualified Data.Array.IArray as A-import Control.Applicative -import Biobase.SubstMatrix+import           Control.Applicative+import           Data.Char (toLower)+import qualified Data.Map as M +import           Biobase.Primary.AA (charAA)+import           Biobase.Types.Odds+import           Data.PrimitiveArray hiding (map)+import qualified Biobase.Primary.AA as AA +import           Biobase.SubstMatrix -fromFile :: FilePath -> IO SubstMatrix+++fromFile :: FilePath -> IO (AASubstMat t) fromFile fname = do   (x:xs) <- dropWhile (("#"==).take 1) . lines <$> readFile fname   let cs = map head . words $ x -- should give us the characters encoding an amino acid-  let ss = map (map read . drop 1 . words) $ xs-  return $ A.accumArray (\_ z -> z) (-10) (('*','*'),('Z','Z')) [ ((k1,k2),z) | (k1,s) <- zip cs ss, (k2,z) <- zip cs s ]+  let ss = map (map DLO . map read . drop 1 . words) $ xs+  let xs = [ ((Z:.charAA k1:.charAA k2),z)+           | (k1,s) <- zip cs ss+           , (k2,z) <- zip cs s+           ]+  return . AASubstMat $ fromAssocs (Z:.AA.Stop:.AA.Stop) (Z:.AA.Z:.AA.Z) (DLO $ -999) xs
BiobaseBlast.cabal view
@@ -1,50 +1,111 @@ name:           BiobaseBlast-version:        0.0.0.1-author:         Christian Hoener zu Siederdissen-maintainer:     choener@tbi.univie.ac.at-homepage:       http://www.tbi.univie.ac.at/~choener/-copyright:      Christian Hoener zu Siederdissen, 2013+version:        0.2.0.0+author:         Christian Hoener zu Siederdissen, Florian Eggenhofer+maintainer:     choener@bioinf.uni-leipzig.de+homepage:       https://github.com/choener/BiobaseBlast+bug-reports:    https://github.com/choener/BiobaseBlast/issues+copyright:      Christian Hoener zu Siederdissen, 2013 - 2017 category:       Bioinformatics-synopsis:       BLAST-related tools license:        GPL-3 license-file:   LICENSE build-type:     Simple stability:      experimental-cabal-version:  >= 1.6.0+cabal-version:  >= 1.10.0+tested-with:    GHC == 7.10.3, GHC == 8.0.1, GHC == 8.0.2+synopsis:       BLAST-related tools description:-                This library contains BLAST-related functionality. For now,-                this library is very limited (and Ketil Malde provides other-                BLAST functionality anyway).+                This library contains BLAST-related functionality:                 .-                We do provide parsers for BLOSUM and PAM matrices.+                - Parser for tabular NCBI BLAST+ output+                - Parsers for BLOSUM and PAM matrices.+                - Specialized substitution functions for (in)complete amino+                  acid / nucleotide triplet substitution.+                - Incomplete nucleotide patterns map one or two nucleotides to+                  an amino acid (need for indel editing in the mitochondria of+                  certain species like /p.polycephalum/).                 .-                The matrices can be found here:+                The matrices are currently not provided but can be found here:                 <ftp://ftp.ncbi.nih.gov/blast/matrices/>    extra-source-files:-  changelog+  changelog.md+  README.md    library-  build-depends:-    base >3 && <5           ,-    array                   ,-    containers              ,-    BiobaseXNA >= 0.7.0.2+  build-depends: base             >= 4.7      && < 5.0+               , aeson            >= 1.0+               , attoparsec       >= 0.13+               , binary           >= 0.7+               , bytestring+               , cereal           >= 0.4+               , containers+               , deepseq          >= 1.3+               , directory+               , vector           >= 0.10+               , vector-th-unbox  >= 0.2+               --+               , BiobaseTypes     == 0.1.2.*+               , BiobaseXNA       == 0.9.3.*+               , PrimitiveArray   == 0.8.0.* +  default-language:+    Haskell2010+  default-extensions: DeriveGeneric+                    , MultiParamTypeClasses+                    , TemplateHaskell+                    , TypeFamilies+                    , TypeOperators   exposed-modules:     Biobase.SubstMatrix     Biobase.SubstMatrix.Import+    Biobase.BLAST+    Biobase.BLAST.Import+    Biobase.BLAST.Types    ghc-options:     -O2   +test-suite properties+  type:+    exitcode-stdio-1.0+  main-is:+    properties.hs+  ghc-options:+    -O2 -threaded -rtsopts -with-rtsopts=-N+  hs-source-dirs:+    tests+  default-language:+    Haskell2010+  default-extensions: BangPatterns+                    , OverloadedStrings+                    , ScopedTypeVariables+                    , TemplateHaskell+  build-depends: base+               , bytestring+               , containers+               , filepath+               , split                      >= 0.2.3+               , tasty                      >= 0.11+               , tasty-quickcheck           >= 0.8+               , tasty-silver               >= 3.1.9+               , tasty-th                   >= 0.1+               --+               , BiobaseBlast+++ source-repository head-  type: git-  location: git://github.com/choener/BiobaseBlast+  type:     git+  location: https://github.com/choener/BiobaseBlast++source-repository this+  type:     git+  location: https://github.com/choener/BiobaseBlast/tree/0.2.0.0+  tag:      0.2.0.0
+ README.md view
@@ -0,0 +1,15 @@+[![Build Status](https://travis-ci.org/choener/BiobaseBlast.svg?branch=master)](https://travis-ci.org/choener/BiobaseBlast)++# BiobaseBlast++This library contains BLAST-related functionality.++++#### Contact++Christian Hoener zu Siederdissen  +Leipzig University, Leipzig, Germany  +choener@bioinf.uni-leipzig.de  +http://www.bioinf.uni-leipzig.de/~choener/  +
− changelog 
@@ -1,2 +0,0 @@-0.0.0.1-		* initial commit
+ changelog.md view
@@ -0,0 +1,16 @@+0.2.0.0+-------++- Added parsing functionality for tabular NCBI BLAST+ output++0.1.0.0+-------++- more efficient tables (using PrimitiveArray)+- updated to newest PrimitiveArray & BiobaseXNA+- travis.yml added++0.0.0.1+-------++- initial commit
+ tests/properties.hs view
@@ -0,0 +1,13 @@++-- |+--+-- TODO some unit tests to check if assumed-known values are in the right+-- table cells++module Main where++++main :: IO ()+main = return ()+