bioinformatics-toolkit 0.3.2 → 0.4.0
raw patch · 10 files changed
+64/−51 lines, 10 files
Files
- bioinformatics-toolkit.cabal +2/−2
- src/Bio/Data/Fasta.hs +3/−1
- src/Bio/GO/GREAT.hs +16/−15
- src/Bio/Motif.hs +1/−1
- src/Bio/Seq.hs +25/−23
- src/Bio/Seq/IO.hs +1/−1
- tests/Tests/ChIPSeq.hs +3/−1
- tests/Tests/GREAT.hs +5/−2
- tests/Tests/Motif.hs +3/−1
- tests/Tests/Seq.hs +5/−4
bioinformatics-toolkit.cabal view
@@ -1,12 +1,12 @@ name: bioinformatics-toolkit-version: 0.3.2+version: 0.4.0 synopsis: A collection of bioinformatics tools description: A collection of bioinformatics tools license: MIT license-file: LICENSE author: Kai Zhang maintainer: kai@kzhang.org-copyright: (c) 2014-2017 Kai Zhang+copyright: (c) 2014-2018 Kai Zhang category: Bio build-type: Simple extra-source-files: README.md
src/Bio/Data/Fasta.hs view
@@ -25,7 +25,9 @@ {-# MINIMAL fromFastaRecord #-} instance BioSeq s a => FastaLike (s a) where- fromFastaRecord (_, xs) = fromBS . B.concat $ xs+ fromFastaRecord (_, xs) = case fromBS (B.concat xs) of+ Left err -> error err+ Right x -> x {-# INLINE fromFastaRecord #-} instance FastaLike Motif where
src/Bio/GO/GREAT.hs view
@@ -14,6 +14,7 @@ import qualified Data.IntervalMap as IM import Data.List (foldl', sortBy) import Data.List.Ordered (nubSort)+import Data.Maybe (fromJust, isNothing) import Bio.Data.Bed import Bio.Utils.Misc (readInt)@@ -35,22 +36,22 @@ -- | Given a gene list and the rule, compute the rulatory domain for each gene getRegulatoryDomains :: AssocRule -> [Gene a] -> [(BED3, a)] getRegulatoryDomains _ [] = error "No gene available for domain assignment!"-getRegulatoryDomains (BasalPlusExtension up dw ext) genes = (extendTail r ext, a) : rs+getRegulatoryDomains (BasalPlusExtension up dw ext) genes = zip+ (loop $ [Nothing] ++ map Just basal ++ [Nothing]) names where- (rs, Just (r,a)) = foldl' f ([], Nothing) $ sortBy (compareBed `on` fst) basal- f (acc, Nothing) (b,x) = (acc, Just (extendHead b ext, x))- f (acc, Just (b', x')) (b,x)- | chrom b' /= chrom b = ( (extendTail b' ext, x') : acc- , Just (extendHead b ext, x) )- | chromEnd b' >= chromStart b = ((b',x') : acc, Just (b,x))- | otherwise = let ext' = min ext $ (chromStart b - chromEnd b') `div` 2- in ((extendTail b' ext', x') : acc, Just (extendHead b ext', x))- extendHead (BED3 chr s e) l | s - l >= 0 = BED3 chr (s-l) e- | otherwise = BED3 chr 0 e- extendTail (BED3 chr s e) l = BED3 chr s (e+l)- basal = flip map genes $ \((chr, tss, str), x) ->- if str then (BED3 chr (tss - up) (tss + dw), x)- else (BED3 chr (tss - dw) (tss + up), x)+ loop (a:b:c:rest) = fn a b c : loop (b:c:rest)+ loop _ = []+ fn left (Just (BED3 chr s e)) right = BED3 chr leftPos rightPos+ where+ leftPos+ | isNothing left || chr /= chrom (fromJust left) = max (s - ext) 0+ | otherwise = min s $ max (s - ext) $ chromEnd $ fromJust left+ rightPos+ | isNothing right || chr /= chrom (fromJust right) = e + ext -- TODO: bound check+ | otherwise = max e $ min (e + ext) $ chromStart $ fromJust right+ (basal, names) = unzip $ sortBy (compareBed `on` fst) $ flip map genes $+ \((chr, tss, str), x) -> if str then (BED3 chr (tss - up) (tss + dw), x)+ else (BED3 chr (tss - dw) (tss + up), x) getRegulatoryDomains _ _ = undefined {-# INLINE getRegulatoryDomains #-}
src/Bio/Motif.hs view
@@ -104,7 +104,7 @@ -- | Convert pwm to consensus sequence, see D. R. Cavener (1987). toIUPAC :: PWM -> DNA IUPAC-toIUPAC (PWM _ pwm) = fromBS . B.pack . map f $ M.toRows pwm+toIUPAC (PWM _ pwm) = unsafeFromBS . B.pack . map f $ M.toRows pwm where f v | snd a > 0.5 && snd a > 2 * snd b = fst a | snd a + snd b > 0.75 = iupac (fst a, fst b)
src/Bio/Seq.hs view
@@ -1,6 +1,7 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables #-} module Bio.Seq ( -- * Alphabet@@ -23,6 +24,7 @@ import Data.Char8 (toUpper) import qualified Data.HashMap.Strict as M import qualified Data.HashSet as S+import Data.Proxy (Proxy (..)) import Prelude hiding (length) -- | Alphabet defined by http://www.chem.qmul.ac.uk/iupac/@@ -54,54 +56,54 @@ class BioSeq' s where toBS :: s a -> B.ByteString+ unsafeFromBS :: B.ByteString -> s a slice :: Int -> Int -> s a -> s a length :: s a -> Int length = B.length . toBS- {-# MINIMAL toBS, slice #-}-+ {-# MINIMAL toBS, slice, unsafeFromBS #-} instance BioSeq' DNA where toBS (DNA s) = s+ unsafeFromBS = DNA slice i l (DNA s) = DNA . B.take l . B.drop i $ s instance BioSeq' RNA where toBS (RNA s) = s+ unsafeFromBS = RNA slice i l (RNA s) = RNA . B.take l . B.drop i $ s instance BioSeq' Peptide where toBS (Peptide s) = s+ unsafeFromBS = Peptide slice i l (Peptide s) = Peptide . B.take l . B.drop i $ s -class BioSeq' s => BioSeq s a where- alphabet :: s a -> S.HashSet Char- fromBS :: B.ByteString -> s a+class BioSeq' seq => BioSeq seq alphabet where+ alphabet :: Proxy (seq alphabet) -> S.HashSet Char+ fromBS :: B.ByteString -> Either String (seq alphabet)+ fromBS input = case B.mapAccumL fun Nothing input of+ (Nothing, r) -> Right $ unsafeFromBS r+ (Just e, _) -> Left $ "Bio.Seq.fromBS: unknown character: " ++ [e]+ where+ fun (Just e) x = (Just e, x)+ fun Nothing x = let x' = toUpper x+ in if x' `S.member` alphabet (Proxy :: Proxy (seq alphabet))+ then (Nothing, x')+ else (Just x', x')+ {-# MINIMAL alphabet #-} instance BioSeq DNA Basic where alphabet _ = S.fromList "ACGT"- fromBS = DNA . B.map (f . toUpper)- where- f x | x `S.member` alphabet (undefined :: DNA Basic) = x- | otherwise = error $ "Bio.Seq.fromBS: unknown character: " ++ [x] instance BioSeq DNA IUPAC where alphabet _ = S.fromList "ACGTNVHDBMKWSYR"- fromBS = DNA . B.map (f . toUpper)- where- f x | x `S.member` alphabet (undefined :: DNA IUPAC) = x- | otherwise = error $ "Bio.Seq.fromBS: unknown character: " ++ [x] instance BioSeq DNA Ext where- alphabet = undefined- fromBS = undefined+ alphabet _ = undefined instance BioSeq RNA Basic where alphabet _ = S.fromList "ACGU"- fromBS = RNA . B.map (f . toUpper)- where- f x | x `S.member` alphabet (undefined :: RNA Basic) = x- | otherwise = error $ "Bio.Seq.fromBS: unknown character: " ++ [x] -- | O(n) Reverse complementary of DNA sequence. rc :: DNA alphabet -> DNA alphabet@@ -112,7 +114,7 @@ 'C' -> 'G' 'G' -> 'C' 'T' -> 'A'- _ -> x+ _ -> x -- | O(n) Compute GC content. gcContent :: DNA alphabet -> Double@@ -130,13 +132,13 @@ 'B' -> x + 0.75 'S' -> x + 1 'W' -> x- _ -> x + 0.5 -- "NMKYR"+ _ -> x + 0.5 -- "NMKYR" in (x', n+1) -- | O(n) Compute single nucleotide frequency.-nucleotideFreq :: BioSeq DNA a => DNA a -> M.HashMap Char Int+nucleotideFreq :: forall a . BioSeq DNA a => DNA a -> M.HashMap Char Int nucleotideFreq dna = B.foldl' f m0 . toBS $ dna where- m0 = M.fromList . zip (S.toList $ alphabet dna) . repeat $ 0+ m0 = M.fromList . zip (S.toList $ alphabet (Proxy :: Proxy (DNA a))) . repeat $ 0 f m x = M.adjust (+1) x m {-# INLINE nucleotideFreq #-}
src/Bio/Seq/IO.hs view
@@ -70,7 +70,7 @@ ">" ++ show chrSize else do hSeek h AbsoluteSeek $ fromIntegral $ headerSize + chrStart + start- (Right . fromBS) <$> B.hGet h (end - start)+ fromBS <$> B.hGet h (end - start) _ -> return $ Left $ "Bio.Seq.getSeq: Cannot find " ++ show chr {-# INLINE getSeq #-}
tests/Tests/ChIPSeq.hs view
@@ -17,9 +17,10 @@ tests :: TestTree tests = testGroup "Test: Bio.ChIPSeq"- [ testCase "rpkm" testRPKM+ [ ] +{- testRPKM :: Assertion testRPKM = do regions <- peaks r1 <- tags $$ rpkmBed regions@@ -27,3 +28,4 @@ let r1' = map (printf "%0.6f") . V.toList $ r1 :: [String] r2' = map (printf "%0.6f") r2 r1' @=? r2'+ -}
tests/Tests/GREAT.hs view
@@ -16,12 +16,14 @@ , (("chr1", 6000, True), 1) , (("chr2", 6000, True), 2) , (("chr2", 2000000, False), 3)+ , (("chr2", 2020000, True), 4) ] results = [ (BED3 "chr1" 0 2000, 0) , (BED3 "chr1" 1000 1007000, 1)- , (BED3 "chr2" 0 1003000, 2)- , (BED3 "chr2" 1003000 3005000, 3)+ , (BED3 "chr2" 0 1007000, 2)+ , (BED3 "chr2" 999000 2015000, 3)+ , (BED3 "chr2" 2005000 3021000, 4) ] loops = [ (BED3 "chr1" 900 1200, BED3 "chr1" 10000 11000)@@ -35,6 +37,7 @@ , (BED3 "chr1" 10000 11000, 0) , (BED3 "chr1" 10000 11000, 1) , (BED3 "chr2" 0 100, 2)+ , (BED3 "chr2" 2015000 2021000,4) ] tests :: TestTree
tests/Tests/Motif.hs view
@@ -18,7 +18,9 @@ import qualified Data.Vector.Unboxed as U dna :: DNA Basic-dna = fromBS $ B.pack $ map f $ take 5000 $ randomRs (0, 3) (mkStdGen 2)+dna = case fromBS (B.pack $ map f $ take 5000 $ randomRs (0, 3) (mkStdGen 2)) of+ Left msg -> error msg+ Right x -> x where f :: Int -> Char f x = case x of
tests/Tests/Seq.hs view
@@ -1,10 +1,11 @@ {-# LANGUAGE OverloadedStrings #-} module Tests.Seq (tests) where -import Bio.Seq+import Bio.Seq+import Data.Either import qualified Data.HashMap.Strict as M-import Test.Tasty-import Test.Tasty.HUnit+import Test.Tasty+import Test.Tasty.HUnit tests :: TestTree tests = testGroup "Test: Bio.Seq"@@ -13,7 +14,7 @@ testNuclFreq :: Assertion testNuclFreq = do- let dna = fromBS "ACTTCCCGGGD" :: DNA IUPAC+ let dna = fromRight undefined $ fromBS "ACTTCCCGGGD" :: DNA IUPAC test = M.lookupDefault undefined 'C' $ nucleotideFreq dna expected = 4 test @=? expected