bioinformatics-toolkit 0.1.0 → 0.1.1
raw patch · 7 files changed
+75/−17 lines, 7 files
Files
- bioinformatics-toolkit.cabal +1/−1
- src/Bio/Data/Bed.hs +29/−10
- src/Bio/Motif.hs +1/−1
- src/Bio/Motif/Alignment.hs +3/−0
- src/Bio/Utils/Misc.hs +15/−3
- tests/Tests/Bed.hs +14/−0
- tests/Tests/Motif.hs +12/−2
bioinformatics-toolkit.cabal view
@@ -2,7 +2,7 @@ -- further documentation, see http://haskell.org/cabal/users-guide/ name: bioinformatics-toolkit-version: 0.1.0+version: 0.1.1 synopsis: A collection of bioinformatics tools description: A collection of bioinformatics tools license: MIT
src/Bio/Data/Bed.hs view
@@ -21,6 +21,7 @@ , splitBed , splitBedBySize , splitBedBySizeLeft+ , splitBedBySizeOverlap , Sorted(..) , sortBed , intersectBed@@ -39,6 +40,7 @@ , fetchSeq , fetchSeq' , compareBed+ , convert ) where import Control.Arrow ((***))@@ -59,7 +61,8 @@ import Bio.Seq import Bio.Seq.IO-import Bio.Utils.Misc (binBySizeLeft, binBySize, bins, readInt, readDouble)+import Bio.Utils.Misc ( binBySizeLeft, binBySize, binBySizeOverlap, bins+ , readInt, readDouble ) -- | a class representing BED-like data, e.g., BED3, BED6 and BED12. BED format -- uses 0-based index (see documentation).@@ -195,8 +198,20 @@ e = chromEnd bed {-# INLINE splitBedBySizeLeft #-} +splitBedBySizeOverlap :: BEDLike b+ => Int -- ^ bin size+ -> Int -- ^ overlap size+ -> b -> [b]+splitBedBySizeOverlap k o bed = map (uncurry (asBed chr)) .+ binBySizeOverlap k o $ (s, e)+ where+ chr = chrom bed+ s = chromStart bed+ e = chromEnd bed+{-# INLINE splitBedBySizeOverlap #-}+ -- | a type to imply that underlying data structure is sorted-newtype Sorted b = Sorted {fromSorted :: b}+newtype Sorted b = Sorted {fromSorted :: b} deriving (Show, Read, Eq) compareBed :: (BEDLike b1, BEDLike b2) => b1 -> b2 -> Ordering compareBed x y = compare x' y'@@ -309,7 +324,7 @@ , _name :: !(Maybe B.ByteString) , _score :: !(Maybe Double) , _strand :: !(Maybe Bool) -- ^ True: "+", False: "-"- } deriving (Eq, Show)+ } deriving (Eq, Show, Read) instance Default BED where def = BED@@ -373,17 +388,17 @@ chromEnd = _chromEnd -- | retreive sequences-fetchSeq :: BioSeq DNA a => Genome -> Conduit BED IO (Either String (DNA a))-fetchSeq g = do gH <- lift $ gHOpen g- table <- lift $ readIndex gH+fetchSeq :: (BioSeq DNA a, MonadIO m) => Genome -> Conduit BED m (Either String (DNA a))+fetchSeq g = do gH <- liftIO $ gHOpen g+ table <- liftIO $ readIndex gH conduitWith gH table- lift $ gHClose gH+ liftIO $ gHClose gH where conduitWith h index' = do bed <- await case bed of Just (BED chr start end _ _ isForward) -> do- dna <- lift $ getSeq h index' (chr, start, end)+ dna <- liftIO $ getSeq h index' (chr, start, end) case isForward of Just False -> yield $ fmap rc dna _ -> yield dna@@ -391,13 +406,13 @@ _ -> return () {-# INLINE fetchSeq #-} -fetchSeq' :: BioSeq DNA a => Genome -> [BED] -> IO [Either String (DNA a)]+fetchSeq' :: (BioSeq DNA a, MonadIO m) => Genome -> [BED] -> m [Either String (DNA a)] fetchSeq' g beds = CL.sourceList beds $= fetchSeq g $$ CL.consume {-# INLINE fetchSeq' #-} -- * BED3 format -data BED3 = BED3 !B.ByteString !Int !Int deriving (Eq, Show)+data BED3 = BED3 !B.ByteString !Int !Int deriving (Eq, Show, Read) instance Default BED3 where def = BED3 "" 0 0@@ -416,3 +431,7 @@ chrom (BED3 f1 _ _) = f1 chromStart (BED3 _ f2 _) = f2 chromEnd (BED3 _ _ f3) = f3++convert :: (BEDLike b1, BEDLike b2) => b1 -> b2+convert b = asBed (chrom b) (chromStart b) (chromEnd b)+{-# INLINE convert #-}
src/Bio/Motif.hs view
@@ -333,7 +333,7 @@ let x = "MOTIF " `B.append` nm y = B.pack $ printf "letter-probability matrix: alength= 4 w= %d nsites= %d E= 0" (size pwm) sites z = B.unlines . map (B.unwords . ("":) . map toShortest) . M.toLists . _mat $ pwm- sites | isNothing (_nSites pwm) = 0+ sites | isNothing (_nSites pwm) = 1 | otherwise = fromJust $ _nSites pwm in B.unlines [x,y,z] {-# INLINE toMEME #-}
src/Bio/Motif/Alignment.hs view
@@ -5,6 +5,9 @@ module Bio.Motif.Alignment ( alignment , alignmentBy+ , linPenal+ , quadPenal+ , cubePenal , mergePWM , buildTree , progressiveMerging
src/Bio/Utils/Misc.hs view
@@ -4,6 +4,7 @@ , bins , binBySize , binBySizeLeft+ , binBySizeOverlap ) where import Data.ByteString.Char8 (ByteString)@@ -30,10 +31,21 @@ in zip xs . tail $ xs {-# INLINE binBySize #-} --- | Including leftovers, the last bin may not have desired size.+binBySizeOverlap :: Int -> Int -> (Int, Int) -> [(Int, Int)]+binBySizeOverlap step overlap (start, end)+ | overlap >= step = error "binBySizeOverlap: overlap > step"+ | otherwise = go start+ where+ go i | i + overlap < end = (i, i + step) : go (i + step - overlap)+ | otherwise = []+{-# INLINE binBySizeOverlap #-}++-- | Including leftovers, the last bin will be extended to match desirable size binBySizeLeft :: Int -> (Int, Int) -> [(Int, Int)]-binBySizeLeft step (start, end) = let xs = [start, start + step .. end-1] ++ [end]- in zip xs . tail $ xs+binBySizeLeft step (start, end) = go start+ where+ go i | i < end = (i, i + step) : go (i + step)+ | otherwise = [] {-# INLINE binBySizeLeft #-} -- | divide a given region into k equal size sub-regions, discarding leftovers
tests/Tests/Bed.hs view
@@ -10,6 +10,7 @@ tests :: TestTree tests = testGroup "Test: Bio.Data.Bed" [ testCase "sortBed" sortBedTest+ , testCase "split" splitBedTest ] sortBedTest :: Assertion@@ -17,3 +18,16 @@ let (Sorted actual) = sortBed beds expect <- fmap V.fromList $ readBed' "tests/data/peaks.sorted.bed" expect @=? actual++splitBedTest :: Assertion+splitBedTest = (s1', s2', s3') @=? (s1, s2, s3)+ where+ bed = BED3 "chr1" 0 99+ s1 = splitBedBySize 20 bed+ s1' = map f [(0, 20), (20, 40), (40, 60), (60, 80)]+ s2 = splitBedBySizeLeft 20 bed+ s2' = map f [(0, 20), (20, 40), (40, 60), (60, 80), (80, 100)]+ s3 = splitBedBySizeOverlap 20 10 bed+ s3' = map f [ ( 0, 20), (10, 30), (20, 40), (30, 50), (40, 60)+ , (50, 70), (60, 80), (70, 90), (80, 100) ]+ f (a,b) = asBed "chr1" a b
tests/Tests/Motif.hs view
@@ -33,6 +33,7 @@ tests = testGroup "Test: Bio.Motif" [ --testCase "IUPAC converting" toIUPACTest testCase "TFBS scanning" findTFBSTest+ , testCase "Max matching score" maxScTest , testCase "pValue calculation" pValueTest ] @@ -53,9 +54,18 @@ actual <- findTFBSSlow def pwm dna (0.6 * optimalScore def pwm) $$ CL.consume assertEqual "findTFBS" expect actual +maxScTest :: Assertion+maxScTest = do+ ms <- motifs+ let expect = map (\x -> maximum $ scores def (_pwm x) dna) ms+ actual = map (\x -> maxMatchSc def (_pwm x) dna) ms+ assertEqual "maxMatchSc" expect actual+ pValueTest :: Assertion pValueTest = do ms <- motifs- let expect = map (pValueToScoreExact 1e-4 def . _pwm) ms- actual = map (pValueToScore 1e-4 def . _pwm) ms+ let expect = map (approx . pValueToScoreExact 1e-4 def . _pwm) ms+ actual = map (approx . pValueToScore 1e-4 def . _pwm) ms assertEqual "pValueToScore" expect actual+ where+ approx x = round $ 10 * x