packages feed

bioinformatics-toolkit 0.5.0 → 0.5.1

raw patch · 9 files changed

+299/−43 lines, 9 filesdep ~base

Dependency ranges changed: base

Files

README.md view
@@ -1,2 +1,2 @@-the-bioinformatician-s-toolkit-==============================+Bioinformatics Algorithms+=========================
bioinformatics-toolkit.cabal view
@@ -1,5 +1,5 @@ name:                bioinformatics-toolkit-version:             0.5.0+version:             0.5.1 synopsis:            A collection of bioinformatics tools description:         A collection of bioinformatics tools license:             MIT@@ -31,6 +31,7 @@     Bio.Data.Bed.Types     Bio.Data.Bam     Bio.Data.Fasta+    Bio.Data.Fastq     Bio.GO     Bio.GO.Parser     Bio.Motif@@ -43,6 +44,7 @@     Bio.RealWorld.GENCODE     Bio.RealWorld.ID     Bio.RealWorld.UCSC+    Bio.RealWorld.Uniprot     Bio.Seq     Bio.Seq.IO     Bio.Utils.Functions@@ -51,7 +53,7 @@     Bio.Utils.Types    build-depends:-      base >=4.8 && <5.0+      base >=4.11 && <5.0     , aeson     , aeson-pretty     , bytestring >=0.10
src/Bio/Data/Bed/Types.hs view
@@ -1,7 +1,25 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TemplateHaskell #-} -module Bio.Data.Bed.Types where+module Bio.Data.Bed.Types+    ( BEDLike(..)+    , BEDConvert(..)+    , BED(..)+    , BED3(..)+    , NarrowPeak(..)+    , npSignal+    , npPvalue+    , npQvalue+    , npPeak+    , BroadPeak(..)+    , bpSignal+    , bpPvalue+    , bpQvalue+    , BEDExt(..)+    , _bed+    , _data+    , BEDTree+    ) where  import           Control.Lens import qualified Data.ByteString.Char8             as B@@ -15,6 +33,20 @@  import           Bio.Utils.Misc                    (readDouble, readInt) +readDoubleNonnegative :: B.ByteString -> Maybe Double+readDoubleNonnegative x | v < 0 = Nothing+                        | otherwise = Just v+  where+    v = readDouble x+{-# INLINE readDoubleNonnegative #-}++readIntNonnegative :: B.ByteString -> Maybe Int+readIntNonnegative x | v < 0 = Nothing+                     | otherwise = Just v+  where+    v = readInt x+{-# INLINE readIntNonnegative #-}+ -- | A class representing BED-like data, e.g., BED3, BED6 and BED12. BED format -- uses 0-based index (see documentation). class BEDLike b where@@ -177,9 +209,9 @@         (readDouble e)         (if f == "." then Nothing else if f == "+" then Just True else Just False)         (readDouble g)-        (if readDouble h < 0 then Nothing else Just $ readDouble h)-        (if readDouble i < 0 then Nothing else Just $ readDouble i)-        (if readInt j < 0 then Nothing else Just $ readInt j)+        (readDoubleNonnegative h)+        (readDoubleNonnegative i)+        (readIntNonnegative j)       where         (a:b:c:d:e:f:g:h:i:j:_) = B.split '\t' l     {-# INLINE fromLine #-}@@ -199,6 +231,62 @@      convert bed = NarrowPeak (bed^.chrom) (bed^.chromStart) (bed^.chromEnd) (bed^.name)         (fromMaybe 0 $ bed^.score) (bed^.strand) 0 Nothing Nothing Nothing++-- | ENCODE broadPeak format: https://genome.ucsc.edu/FAQ/FAQformat.html#format13+data BroadPeak = BroadPeak+    { _bpChrom  :: !B.ByteString+    , _bpStart  :: !Int+    , _bpEnd    :: !Int+    , _bpName   :: !(Maybe B.ByteString)+    , _bpScore  :: !Double+    , _bpStrand :: !(Maybe Bool)+    , _bpSignal  :: !Double+    , _bpPvalue :: !(Maybe Double)+    , _bpQvalue :: !(Maybe Double)+    } deriving (Eq, Show, Read)++makeLensesFor [ ("_bpSignal", "bpSignal")+              , ("_bpPvalue", "bpPvalue")+              , ("_bpQvalue", "bpQvalue")+              ] ''BroadPeak++instance BEDLike BroadPeak where+    chrom = lens _bpChrom (\bed x -> bed { _bpChrom = x })+    chromStart = lens _bpStart (\bed x -> bed { _bpStart = x })+    chromEnd = lens _bpEnd (\bed x -> bed { _bpEnd = x })+    name = lens _bpName (\bed x -> bed { _bpName = x })+    score = lens (Just . _bpScore) (\bed x -> bed { _bpScore = fromJust x })+    strand = lens _bpStrand (\bed x -> bed { _bpStrand = x })++instance BEDConvert BroadPeak where+    asBed chr s e = BroadPeak chr s e Nothing 0 Nothing 0 Nothing Nothing++    fromLine l = BroadPeak a (readInt b) (readInt c)+        (if d == "." then Nothing else Just d)+        (readDouble e)+        (if f == "." then Nothing else if f == "+" then Just True else Just False)+        (readDouble g)+        (readDoubleNonnegative h)+        (readDoubleNonnegative i)+      where+        (a:b:c:d:e:f:g:h:i:_) = B.split '\t' l+    {-# INLINE fromLine #-}++    toLine (BroadPeak a b c d e f g h i) = B.intercalate "\t"+        [ a, fromJust $ packDecimal b, fromJust $ packDecimal c, fromMaybe "." d+        , toShortest e+        , case f of+            Nothing   -> "."+            Just True -> "+"+            _         -> "-"+        , toShortest g, fromMaybe "-1" $ fmap toShortest h+        , fromMaybe "-1" $ fmap toShortest i+        ]+    {-# INLINE toLine #-}++    convert bed = BroadPeak (bed^.chrom) (bed^.chromStart) (bed^.chromEnd) (bed^.name)+        (fromMaybe 0 $ bed^.score) (bed^.strand) 0 Nothing Nothing+  data BEDExt bed a = BEDExt     { _ext_bed :: bed
+ src/Bio/Data/Fastq.hs view
@@ -0,0 +1,107 @@+module Bio.Data.Fastq+    ( Fastq(..)+    , parseFastqC+    , parseFastqUnsafeC+    , fastqToByteString+    , trimPolyA+    ) where++import           Conduit+import           Control.Monad         (when)+import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString as BS+import           Data.Maybe            (isJust)++-- | A FASTQ file normally uses four lines per sequence.+--+--     * Line 1 begins with a '@' character and is followed by a sequence+--       identifier and an optional description (like a FASTA title line).+--+--     * Line 2 is the raw sequence letters.+--+--     * Line 3 begins with a '+' character and is optionally followed by the+--       same sequence identifier (and any description) again.+--+--     * Line 4 encodes the quality values for the sequence in Line 2, and must+--       contain the same number of symbols as letters in the sequence.+data Fastq = Fastq+    { fastqSeqId   :: B.ByteString+    , fastqSeq     :: B.ByteString+    , fastqSeqInfo :: B.ByteString+    , fastqSeqQual :: B.ByteString+    } deriving (Show, Eq)++parseFastqC :: Monad m => ConduitT B.ByteString Fastq m ()+parseFastqC = linesUnboundedAsciiC .| conduit+  where+    conduit = do+        l1 <- await+        l2 <- await+        l3 <- await+        l4 <- await+        case mkFastqRecord <$> l1 <*> l2 <*> l3 <*> l4 of+            Nothing -> when (isJust l1) $ error "file ends prematurely"+            Just x  -> yield x >> conduit+{-# INLINE parseFastqC #-}++parseFastqUnsafeC :: Monad m => ConduitT B.ByteString Fastq m ()+parseFastqUnsafeC = linesUnboundedAsciiC .| conduit+  where+    conduit = do+        l1 <- await+        l2 <- await+        l3 <- await+        l4 <- await+        case mkFastqRecordUnsafe <$> l1 <*> l2 <*> l3 <*> l4 of+            Nothing -> when (isJust l1) $ error "file ends prematurely"+            Just x  -> yield x >> conduit+{-# INLINE parseFastqUnsafeC #-}++fastqToByteString :: Fastq -> [B.ByteString]+fastqToByteString (Fastq a b c d) = ['@' `B.cons` a, b, '+' `B.cons` c, d]+{-# INLINE fastqToByteString #-}++-- | Make Fastq record from Bytestrings, without sanity check.+mkFastqRecordUnsafe :: B.ByteString   -- ^ First line+                    -> B.ByteString   -- ^ Second line+                    -> B.ByteString   -- ^ Third line+                    -> B.ByteString   -- ^ Fourth line+                    -> Fastq+mkFastqRecordUnsafe l1 l2 l3 l4 = Fastq (B.tail l1) l2 (B.tail l3) l4+{-# INLINE mkFastqRecordUnsafe #-}++mkFastqRecord :: B.ByteString   -- ^ First line+              -> B.ByteString   -- ^ Second line+              -> B.ByteString   -- ^ Third line+              -> B.ByteString   -- ^ Fourth line+              -> Fastq+mkFastqRecord l1 l2 l3 l4 = Fastq (parseLine1 l1) (parseLine2 l2)+    (parseLine3 l3) (parseLine4 l4)+  where+    parseLine1 x | B.head x == '@' = B.tail x+                 | otherwise = error $ "Parse line 1 failed: " ++ B.unpack x+    parseLine2 x | B.all f x = x+                 | otherwise = error $ "Parse line 2 failed: " ++ B.unpack x+      where+        f 'C' = True+        f 'G' = True+        f 'T' = True+        f 'A' = True+        f 'N' = True+        f _   = False+    parseLine3 x | B.head x == '+' = B.tail x+                 | otherwise = error $ "Parse line 3 failed: " ++ B.unpack x+    parseLine4 x | BS.all f x = x+                 | otherwise = error $ "Parse line 4 failed: " ++ B.unpack x+      where+        f b = let b' = fromIntegral b :: Int+              in b' >= 33 && b' <= 126+{-# INLINE mkFastqRecord #-}++-- | Remove trailing 'A'+trimPolyA :: Int -> Fastq -> Fastq+trimPolyA n f@(Fastq a b c d)+    | B.length trailing >= n = Fastq a b' c $ B.take (B.length b') d+    | otherwise = f+  where+    (b', trailing) = B.spanEnd (=='A') b
src/Bio/RealWorld/ENCODE.hs view
@@ -35,6 +35,7 @@ import qualified Data.Sequence as S import qualified Data.Text as T import qualified Data.Vector as V+import Data.Semigroup (Semigroup(..)) import Network.HTTP.Conduit import Data.Default.Class @@ -54,9 +55,8 @@         g y' | S.null y' = ""              | otherwise =  foldr1 (\a b -> b ++ ('&':a)) y' -instance Monoid KeyWords where-    mempty = KeyWords S.empty S.empty-    mappend (KeyWords a b) (KeyWords a' b') = KeyWords (a S.>< a') (b S.>< b')+instance Semigroup KeyWords where+    (<>) (KeyWords a b) (KeyWords a' b') = KeyWords (a S.>< a') (b S.>< b')  base :: String base = "https://www.encodeproject.org/"
+ src/Bio/RealWorld/Uniprot.hs view
@@ -0,0 +1,36 @@+{-# LANGUAGE OverloadedStrings #-}++module Bio.RealWorld.Uniprot+    ( mapID+    ) where++import           Conduit+import qualified Data.ByteString.Char8 as B+import qualified Data.HashMap.Strict   as M+import           Network.HTTP.Conduit++base :: String+base = "http://www.uniprot.org/uploadlists/"++mapID :: [B.ByteString]   -- ^ A list of IDs+      -> B.ByteString     -- ^ From database+      -> B.ByteString     -- ^ To database+      -> IO [Maybe B.ByteString]+mapID ids from to = do+    initReq <- parseRequest base+    let request = setQueryString query initReq+            { method = "GET"+            , requestHeaders = [("User-Agent", "kk@test.org")]+            }+    manager <- newManager tlsManagerSettings+    r <- fmap M.fromList $ runResourceT $ do+        response <- http request manager+        runConduit $ responseBody response .| linesUnboundedAsciiC .|+            (dropC 1 >> mapC ((\[a,b] -> (a,b)) . B.split '\t')) .| sinkList+    return $ map (flip M.lookup r) ids+  where+    query = [ ("from", Just from)+            , ("to", Just to)+            , ("format", Just "tab")+            , ("query", Just $ B.unwords ids)+            ]
src/Bio/Seq.hs view
@@ -49,9 +49,11 @@ instance Show (DNA a) where     show (DNA s) = B.unpack s +instance Semigroup (DNA a) where+    (<>) (DNA x) (DNA y) = DNA (x <> y)+ instance Monoid (DNA a) where     mempty = DNA B.empty-    mappend (DNA x) (DNA y) = DNA (x `mappend` y)     mconcat dnas = DNA . B.concat . map toBS $ dnas  class BioSeq' s where
src/Bio/Utils/Functions.hs view
@@ -17,12 +17,13 @@ ) where  import Data.Bits (shiftR)-import Data.List (foldl')+import Data.List (foldl', groupBy)+import Data.Function (on) import Data.Ord (comparing) import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U-import qualified Data.Matrix.Unboxed as MU+import qualified Data.Matrix as M import Statistics.Sample (meanVarianceUnb, mean) import Statistics.Function (sortBy) @@ -155,19 +156,18 @@ -}  -- | Columns are samples, rows are features / genes.--- TODO: handle ties.-quantileNormalization :: MU.Matrix Double -> MU.Matrix Double-quantileNormalization mat = fromColumns $-    map (snd . (U.unzip :: U.Vector (Int, Double) -> (U.Vector Int, U.Vector Double)) . sortBy (comparing fst)) $ MU.toColumns $-    fromRows $ map f $ MU.toRows $ fromColumns $-    map (sortBy (comparing snd) . U.zip (U.enumFromN 0 n)) $ MU.toColumns mat+quantileNormalization :: M.Matrix Double -> M.Matrix Double+quantileNormalization mat = M.fromColumns $ map+    (fst . G.unzip . sortBy (comparing snd) . G.fromList . concatMap f .+    groupBy ((==) `on` (snd . snd)) . zip averages . G.toList) $+    M.toColumns srtMat   where-    n = MU.rows mat-    f xs = let m = mean $ snd $ U.unzip xs-           in U.map (\(i,_) -> (i,m)) xs--fromRows :: U.Unbox a => [U.Vector a] -> MU.Matrix a-fromRows = MU.fromRows--fromColumns :: U.Unbox a => [U.Vector a] -> MU.Matrix a-fromColumns = MU.fromColumns+    f [(a,(b,c))] = [(a,b)]+    f xs = let m = mean $ U.fromList $ fst $ unzip xs+           in map (\(_,(i,_)) -> (m, i)) xs+    srtMat :: M.Matrix (Int, Double)+    srtMat = M.fromColumns $ map (sortBy (comparing snd) . G.zip (G.enumFromN 0 n)) $+        M.toColumns mat+    averages = map (mean . snd . G.unzip) $ M.toRows srtMat+    n = M.rows mat+{-# INLINE quantileNormalization #-}
tests/Tests/Tools.hs view
@@ -3,23 +3,44 @@ import Bio.Utils.Functions import Test.Tasty import Test.Tasty.HUnit-import qualified Data.Matrix.Unboxed as MU+import qualified Data.Matrix as MU+import Text.Printf (printf)  tests :: TestTree tests = testGroup "Test: Bio.Utils"-    [ testCase "quantileNormalization" quantileNormalizationTest+    [ quantileNormalizationTest     ] -quantileNormalizationTest :: Assertion-quantileNormalizationTest = after @=? quantileNormalization before+quantileNormalizationTest :: TestTree+quantileNormalizationTest = testGroup "quantile normalization"+    [ testCase "case 1" $ y1 @=? quantileNormalization x1+    , testCase "case 2" $ MU.map (printf "%.2f") y2 @=?+        (MU.map (printf "%.2f") (quantileNormalization x2) :: MU.Matrix String)+    ]   where-    before = MU.fromLists [ [2, 4, 4, 5]-                          , [5, 14, 4, 7]-                          , [4, 8, 6, 9]-                          , [3, 8, 5, 8]-                          , [3, 9, 3, 5] ]-    after = MU.fromLists [ [3.5, 3.5, 5.0, 5.0]-                         , [8.5, 8.5, 5.5, 5.5]-                         , [6.5, 5.0, 8.5, 8.5]-                         , [5.0, 5.5, 6.5, 6.5]-                         , [5.5, 6.5, 3.5, 3.5] ]+    x1 :: MU.Matrix Double+    x1 = MU.fromLists+        [ [2,  4, 4, 5]+        , [5, 14, 4, 7]+        , [4,  8, 6, 9]+        , [3,  8, 5, 8]+        , [3,  9, 3, 5] ]+    y1 :: MU.Matrix Double+    y1 = MU.fromLists+        [ [ 3.5,  3.5, 5.25, 4.25]+        , [ 8.5,  8.5, 5.25,  5.5]+        , [ 6.5, 5.25,  8.5,  8.5]+        , [5.25, 5.25,  6.5,  6.5]+        , [5.25,  6.5,  3.5, 4.25] ]+    x2 :: MU.Matrix Double+    x2 = MU.fromLists+        [ [5, 4, 3]+        , [2, 1, 4]+        , [3, 4, 6]+        , [4, 2, 8] ]+    y2 :: MU.Matrix Double+    y2 = MU.fromLists+        [ [ 5.666667,  5.166667, 2 ]+        , [ 2,  2, 3]+        , [3, 5.166667,  4.666667]+        , [4.666667,  3,  5.666667] ]