diff --git a/app/mergeMotifs.hs b/app/mergeMotifs.hs deleted file mode 100644 --- a/app/mergeMotifs.hs +++ /dev/null @@ -1,229 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} - -import AI.Clustering.Hierarchical hiding (drawDendrogram) -import Control.Monad (forM, forM_) -import qualified Data.ByteString.Char8 as B -import Data.Default.Class -import Data.Double.Conversion.ByteString (toShortest) -import Data.List -import Data.List.Split (splitOn) -import Data.Ord -{- -import Diagrams.Backend.Cairo -import Diagrams.Plots.Dendrogram -import Diagrams.Prelude (dims2D, strutX, (|||)) --} -import Options.Applicative -import System.IO -import Text.Printf - -import Bio.Data.Fasta -import Bio.Motif -import Bio.Motif.Alignment -import Bio.Motif.Merge -import Bio.Seq (toBS) -import Bio.Utils.Functions - -data CMD = Merge { mergeInput :: FilePath - , mode :: String - , thres :: Double - , alignOption :: AlignOption - , mergeOutput :: FilePath - } - | Dist { inputA :: FilePath - , inputB :: Maybe FilePath - , alignOption :: AlignOption - } - | Cluster { clusterInput :: !FilePath - , cutoff :: !Double - , alignOption :: !AlignOption - } - -mergeParser :: Parser CMD -mergeParser = Merge - <$> argument str (metavar "INPUT") - <*> strOption - ( long "mode" - <> short 'm' - <> value "iter" - <> metavar "MODE" - <> help "Merging algorithm, could be iter or tree, default is iter" ) - <*> option auto - ( long "thres" - <> short 't' - <> value 0.2 - <> metavar "THRESHOLD" - <> help "two motifs that have distance belowing threshold would be merged, default is 0.2" ) - <*> alignParser - <*> strOption - ( long "output" - <> short 'o' - <> value "merged_output.meme" - <> metavar "OUTPUT" ) - -distParser :: Parser CMD -distParser = Dist - <$> strOption - ( short 'a' - <> metavar "INPUT_A" ) - <*> (optional . strOption) - ( short 'b' - <> metavar "INPUT_B" ) - <*> alignParser - -clusterParser :: Parser CMD -clusterParser = Cluster - <$> argument str (metavar "INPUT") - <*> option auto - ( long "height" - <> short 'h' - <> value 0.2 - <> metavar "HEIGHT" - <> help "Cut hierarchical tree at given height. Default: 0.2" ) - <*> alignParser - -data AlignOption = AlignOption - { gap :: Double - , gapMode :: String - , avgMode :: String - } - -alignParser :: Parser AlignOption -alignParser = AlignOption - <$> option auto - ( long "gap" - <> short 'g' - <> value 0.05 - <> metavar "GAP_PENALTY" - <> help "Gap penalty, default: 0.05" ) - <*> strOption - ( long "gap_mode" - <> value "exp" - <> metavar "GAP_MODE" - <> help "Gap penalty mode, one of linear, quadratic, cubic, and exp. default: exp." ) - <*> strOption - ( long "avg_mode" - <> value "l1" - <> metavar "AVERAGE_MODE" - <> help "Averaging function, one of l1, l2, l3, max. default: l1." ) - -treeMerge :: Double -> String -> [Motif] -> AlignOption - -> ([Motif], Dendrogram Motif) -treeMerge th pre ms alignOpt = (zipWith f [0::Int ..] $ map merge $ tree `cutAt` th, tree) - where - f i (suffix, pwm) = Motif ((B.pack $ pre ++ "_" ++ show i ++ "_" ++ show (toIUPAC pwm)) - `B.append` "(" - `B.append` suffix - `B.append` ")" ) pwm - merge tr = ( B.intercalate "+" $ map _name $ flatten tr - , dilute $ mergeTreeWeighted align tr) - tree = buildTree align ms - align = mkAlignFn alignOpt -{-# INLINE treeMerge #-} - -getSuffix :: String -> String -getSuffix = last . splitOn "." -{-# INLINE getSuffix #-} - -mkAlignFn :: AlignOption -> AlignFn -mkAlignFn (AlignOption gap gapMode avgMode) = alignmentBy jsd (pFn gap) avgFn - where - pFn = case gapMode of - "linear" -> linPenal - "quadratic" -> quadPenal - "cubic" -> cubPenal - "exp" -> expPenal - _ -> error "Unknown gap mode" - avgFn = case avgMode of - "l1" -> l1 - "l2" -> l2 - "l3" -> l3 - "max" -> lInf - _ -> error "Unknown average mode" -{-# INLINE mkAlignFn #-} - -readMotif :: FilePath -> IO [Motif] -readMotif fl = case getSuffix fl of - "fasta" -> readFasta' fl - _ -> readMEME fl - -writeMotif :: FilePath -> [Motif] -> IO () -writeMotif fl motifs = case getSuffix fl of - "fasta" -> writeFasta fl motifs - _ -> writeMEME fl motifs def - - -defaultMain :: CMD -> IO () -defaultMain (Dist a b alignOpt) = do - motifsA <- readMotif a - pairs <- case b of - Nothing -> return $ comb motifsA - Just b' -> do - motifsB <- readMotif b' - return [(x,y) | x <- motifsA, y <- motifsB] - forM_ pairs $ \(x,y) -> do - let (d,_) = alignFn (_pwm x) $ (_pwm y) - B.putStrLn $ B.intercalate "\t" [_name x, _name y, toShortest d] - where - alignFn = mkAlignFn alignOpt - -defaultMain (Merge inFl m th alignOpt outFl)= do - motifs <- readMotif inFl - let motifNumber = length motifs - - hPutStrLn stderr $ printf "Merging Mode: %s" m - hPutStrLn stderr $ printf "Read %d motifs" motifNumber - - motifs' <- case m of - "tree" -> do - let (newMotifs, tree) = treeMerge th "" motifs alignOpt - fn x = B.unpack (_name x) ++ ": " ++ B.unpack (toBS $ toIUPAC $ _pwm x) - - {- - case svg of - Just fl -> do - let w = 80 - h = 5 * fromIntegral motifNumber - renderCairo fl (dims2D (10*w) (10*h)) $ drawDendrogram w h th tree fn ||| strutX 40 - return newMotifs - -} - return newMotifs - "iter" -> do - let rs = iterativeMerge (mkAlignFn alignOpt) th motifs - forM rs $ \(nm, pwm, ws) -> do - let pwm' = dilute (pwm, ws) - return $ Motif (B.intercalate "+" nm) pwm - -- _ -> error "Unkown mode!" - - hPutStrLn stderr $ printf "Write %d motifs" (length motifs') - - writeMotif outFl motifs' - -defaultMain (Cluster inFl c alignOpt) = do - motifs <- readMotif inFl - let tree = buildTree align motifs - align = mkAlignFn alignOpt - forM_ (tree `cutAt` c) $ \t -> - B.putStrLn $ B.intercalate "\t" $ map _name $ flatten t -{-# INLINE defaultMain #-} - -comb :: [a] -> [(a,a)] -comb (y:ys) = zip (repeat y) ys ++ comb ys -comb _ = [] -{-# INLINE comb #-} - -main :: IO () -main = execParser opts >>= defaultMain - where - opts = info (helper <*> parser) - ( fullDesc - <> header (printf "Compare, align and merge motifs, version %s" v)) - v = "0.2.0" :: String - parser = subparser $ ( - command "merge" (info (helper <*> mergeParser) $ - fullDesc <> progDesc "Merge motifs") - <> command "dist" (info (helper <*> distParser) $ - fullDesc <> progDesc "Align and compute pairwise distances") - <> command "cluster" (info (helper <*> clusterParser) $ - fullDesc <> progDesc "Perform hierarchical clustering on motifs") - ) diff --git a/app/mkindex.hs b/app/mkindex.hs deleted file mode 100644 --- a/app/mkindex.hs +++ /dev/null @@ -1,11 +0,0 @@ -module Main where - -import Bio.Seq.IO (mkIndex) -import Shelly -import qualified Data.Text as T -import System.Environment - -main :: IO () -main = do [inDir, outF] <- getArgs - fls <- shelly . lsT . fromText . T.pack $ inDir - mkIndex (map T.unpack fls) outF diff --git a/app/viewSeq.hs b/app/viewSeq.hs deleted file mode 100644 --- a/app/viewSeq.hs +++ /dev/null @@ -1,15 +0,0 @@ -module Main where - -import qualified Data.ByteString.Char8 as B -import Bio.Seq -import Bio.Seq.IO -import System.Environment - -main :: IO () -main = do - [fl, chr, start, end] <- getArgs - withGenome fl $ \g -> do - s <- getSeq g (B.pack chr, read start, read end) :: IO (Either String (DNA IUPAC)) - case s of - Left err -> error err - Right x -> B.putStrLn . toBS $ x diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal --- a/bioinformatics-toolkit.cabal +++ b/bioinformatics-toolkit.cabal @@ -1,8 +1,5 @@ --- Initial bioinformatics-toolkit.cabal generated by cabal init. For --- further documentation, see http://haskell.org/cabal/users-guide/ - name: bioinformatics-toolkit -version: 0.2.4 +version: 0.3.0 synopsis: A collection of bioinformatics tools description: A collection of bioinformatics tools license: MIT @@ -17,7 +14,7 @@ data-files: tests/data/*.bed tests/data/*.sorted.bed --- tests/data/*.bam + tests/data/*.bam -- tests/data/*.bam.bai tests/data/*.meme tests/data/*.fasta @@ -30,8 +27,7 @@ Bio.ChIPSeq Bio.ChIPSeq.FragLen Bio.Data.Bed --- Bio.Data.BigWig --- Bio.Data.Bam + Bio.Data.Bam Bio.Data.Fasta Bio.GO Bio.GO.Parser @@ -53,79 +49,38 @@ Bio.Utils.Overlap Bio.Utils.Types - -- other-modules: - -- other-extensions: build-depends: base >=4.8 && <5.0 , aeson , aeson-pretty --- , bbi - , binary , bytestring >=0.10 , bytestring-lexing >=0.5 , case-insensitive , clustering - , colour , conduit-combinators , containers >=0.5 , data-default-class , double-conversion + , HsHTSLib , http-conduit >=2.1.8 , hexpat + , IntervalMap >=0.5.0.0 , matrices >=0.4.3 , mtl >=2.1.3.1 , math-functions , parallel >=3.2 , primitive - , palette , split , statistics >=0.13.2.1 , text >=0.11 , transformers >=0.3.0.0 , unordered-containers >=0.2 + , word8 , vector , vector-algorithms - , word8 - , IntervalMap >=0.5.0.0 default-language: Haskell2010 -executable mkindex - hs-source-dirs: app - main-is: mkindex.hs - build-depends: - base >=4.8 && <5.0 - , bioinformatics-toolkit - , shelly - , text - default-language: Haskell2010 - -executable viewSeq - hs-source-dirs: app - main-is: viewSeq.hs - build-depends: - base >=4.8 && <5.0 - , bioinformatics-toolkit - , bytestring - default-language: Haskell2010 - -executable mergeMotifs - hs-source-dirs: app - main-is: mergeMotifs.hs - build-depends: - base >=4.6 && <5.0 - , bioinformatics-toolkit - , bytestring - , clustering - , data-default-class - , double-conversion - , split - , optparse-applicative --- , haskell-plot >=0.2.0.0 --- , diagrams-cairo >=1.3 --- , diagrams-lib >=1.3 - default-language: Haskell2010 - benchmark bench type: exitcode-stdio-1.0 main-is: benchmarks/bench.hs @@ -148,6 +103,7 @@ main-is: test.hs other-modules: Tests.Bed + , Tests.Bam , Tests.Motif , Tests.Seq , Tests.GREAT diff --git a/src/Bio/Data/Bam.hs b/src/Bio/Data/Bam.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Data/Bam.hs @@ -0,0 +1,45 @@ +{-# LANGUAGE FlexibleContexts #-} +module Bio.Data.Bam + ( Bam + , HeaderState + , runBam + , readBam + , writeBam + , bamToBed + ) where + +import Bio.Data.Bed +import Bio.HTS +import Bio.HTS.Types (Bam, FileHeader (..)) +import Conduit +import Control.Monad.State (get, lift) + +bamToBed :: Conduit Bam HeaderState BED +bamToBed = mapMC f =$= concatC + where + f bam = do + BamHeader hdr <- lift get + case getChr hdr bam of + Just chr -> + let start = fromIntegral $ position bam + end = fromIntegral $ endPos bam + nm = Just $ qName bam + strand = Just $ not $ isRev bam + in return $ Just $ BED chr start end nm Nothing strand + _ -> return Nothing +{-# INLINE bamToBed #-} + +{- +viewBam :: IdxHandle -> (B.ByteString, Int, Int) -> Source IO Bam1 +viewBam handle (chr, s, e) = case lookupTarget (idxHeader handle) chr of + Nothing -> return () + Just chrId -> do + q <- lift $ query handle chrId (fromIntegral s,fromIntegral e) + go q + where + go q' = do r <- lift $ next q' + case r of + Nothing -> return () + Just bam -> yield bam >> go q' +{-# INLINE viewBam #-} +-} diff --git a/src/Bio/Motif.hs b/src/Bio/Motif.hs --- a/src/Bio/Motif.hs +++ b/src/Bio/Motif.hs @@ -63,31 +63,32 @@ size :: PWM -> Int size (PWM _ mat) = M.rows mat --- | information content of a poistion in pwm +-- | Information content of a poistion in pwm. (Not implemented) ic :: PWM -> Int -> Double ic = undefined --- | extract sub-PWM given starting position and length, zero indexed +-- | Extract sub-PWM given starting position and length, zero indexed. subPWM :: Int -> Int -> PWM -> PWM subPWM i l (PWM n mat) = PWM n $ M.subMatrix (i,0) (i+l,3) mat {-# INLINE subPWM #-} --- | reverse complementary of PWM +-- | Reverse complementary of PWM. rcPWM :: PWM -> PWM rcPWM (PWM n mat) = PWM n . M.fromVector d . U.reverse . M.flatten $ mat where d = M.dim mat {-# INLINE rcPWM #-} --- | GC content of PWM +-- | GC content of PWM. gcContentPWM :: PWM -> Double gcContentPWM (PWM _ mat) = loop 0 0 / fromIntegral m where m = M.rows mat - loop !acc !i | i >= m = acc - | otherwise = - let acc' = acc + M.unsafeIndex mat (i,1) + M.unsafeIndex mat (i,2) - in loop acc' (i+1) + loop !acc !i + | i >= m = acc + | otherwise = + let acc' = acc + M.unsafeIndex mat (i,1) + M.unsafeIndex mat (i,2) + in loop acc' (i+1) data Motif = Motif { _name :: !B.ByteString @@ -95,13 +96,13 @@ } deriving (Show, Read) -- | background model which consists of single nucletide frequencies, and di-nucletide --- frequencies +-- frequencies. newtype Bkgd = BG (Double, Double, Double, Double) instance Default Bkgd where def = BG (0.25, 0.25, 0.25, 0.25) --- | convert pwm to consensus sequence, see D. R. Cavener (1987). +-- | 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 where @@ -121,7 +122,7 @@ sort' (x, y) | x > y = (y, x) | otherwise = (x, y) --- | get scores of a long sequences at each position +-- | Get scores of a long sequences at each position. scores :: Bkgd -> PWM -> DNA a -> [Double] scores bg p@(PWM _ pwm) dna = go $! toBS dna where @@ -130,7 +131,7 @@ len = M.rows pwm {-# INLINE scores #-} --- | a streaming version of scores +-- | A streaming version of scores. scores' :: Monad m => Bkgd -> PWM -> DNA a -> Source m Double scores' bg p@(PWM _ pwm) dna = go 0 where @@ -146,7 +147,7 @@ score bg p dna = scoreHelp bg p $! toBS dna {-# INLINE score #-} --- | the best possible score for a pwm +-- | The best possible matching score of a pwm. optimalScore :: Bkgd -> PWM -> Double optimalScore (BG (a, c, g, t)) (PWM _ pwm) = foldl' (+) 0 . map f . M.toRows $ pwm where f xs = let (i, s) = U.maximumBy (comparing snd) . @@ -190,7 +191,8 @@ pseudoCount = 0.0001 {-# INLINE scoreHelp #-} --- | the probability that a kmer is generated by background (0-orderd model) +-- | The probability that a kmer is generated by background based on a +-- 0-orderd Markov model. pBkgd :: Bkgd -> B.ByteString -> Double pBkgd (BG (a, c, g, t)) = B.foldl' f 1 where @@ -204,14 +206,20 @@ _ -> undefined {-# INLINE pBkgd #-} +-- | calculate the minimum motif mathching score that would produce a kmer with +-- p-Value less than the given number. This score would then be used to search +-- for motif occurrences with significant p-Value +pValueToScore :: Double -> Bkgd -> PWM -> Double +pValueToScore p bg pwm = cdf' (scoreCDF bg pwm) $ 1 - p --- | unlike pValueToScore, this version compute the exact score but much slower --- and is inpractical for long motifs +-- | Unlike pValueToScore, this version compute the exact score but much slower +-- and is inpractical for long motifs. pValueToScoreExact :: Double -- ^ desirable p-Value -> Bkgd -> PWM -> Double -pValueToScoreExact p bg pwm = go 0 0 . sort' . map ((scoreHelp bg pwm &&& pBkgd bg) . B.pack) . replicateM l $ "ACGT" +pValueToScoreExact p bg pwm = go 0 0 $ sort' $ + map ((scoreHelp bg pwm &&& pBkgd bg) . B.pack) $ replicateM l "ACGT" where sort' xs = U.create $ do v <- U.unsafeThaw . U.fromList $ xs I.sortBy (flip (comparing fst)) v @@ -220,13 +228,7 @@ | otherwise = go (acc + snd (vec U.! i)) (i+1) vec l = size pwm --- | calculate the minimum motif mathching score that would produce a kmer with --- p-Value less than the given number. This score would then be used to search --- for motif occurrences with significant p-Value -pValueToScore :: Double -> Bkgd -> PWM -> Double -pValueToScore p bg pwm = cdf' (scoreCDF bg pwm) $ 1 - p - --- | The cumulative distribution function in the form of (x, P(X <= x)) +-- | The cumulative distribution function in the form of (x, P(X <= x)). newtype CDF = CDF (U.Vector (Double, Double)) deriving (Read, Show) -- P(X <= x) @@ -245,7 +247,7 @@ n = U.length v {-# INLINE cdf #-} --- the inverse of cdf +-- | The inverse of cdf. cdf' :: CDF -> Double -> Double cdf' (CDF v) p | p > 1 || p < 0 = error "p must be in [0,1]" @@ -268,7 +270,7 @@ truncateCDF x (CDF v) = CDF $ U.filter ((>=x) . snd) v {-# INLINE truncateCDF #-} --- approximate the cdf of motif matching scores +-- | Approximate the cdf of motif matching scores scoreCDF :: Bkgd -> PWM -> CDF scoreCDF (BG (a,c,g,t)) pwm = toCDF $ loop (U.singleton 1, const 0) 0 where @@ -319,11 +321,11 @@ | otherwise = log x {-# INLINE scoreCDF #-} --- | get pwm from a matrix +-- | Get pwm from a matrix. toPWM :: [B.ByteString] -> PWM toPWM x = PWM Nothing . M.fromLists . map (map readDouble.B.words) $ x --- | pwm to bytestring +-- | Convert pwm to bytestring. fromPWM :: PWM -> B.ByteString fromPWM = B.unlines . map (B.unwords . map toShortest) . M.toLists . _mat diff --git a/src/Bio/Motif/Merge.hs b/src/Bio/Motif/Merge.hs --- a/src/Bio/Motif/Merge.hs +++ b/src/Bio/Motif/Merge.hs @@ -9,6 +9,7 @@ , mergeTreeWeighted , iterativeMerge , buildTree + , cutTreeBy )where import AI.Clustering.Hierarchical hiding (size) diff --git a/src/Bio/RealWorld/BioGRID.hs b/src/Bio/RealWorld/BioGRID.hs --- a/src/Bio/RealWorld/BioGRID.hs +++ b/src/Bio/RealWorld/BioGRID.hs @@ -78,7 +78,7 @@ -- | retreive first 10,000 records fetchByGeneNames :: [String] -> IO [TAB2] fetchByGeneNames genes = do - initReq <- parseUrl $ intercalate "&" [url, geneList, tax, accessKey] + initReq <- parseRequest $ intercalate "&" [url, geneList, tax, accessKey] let request = initReq { method = "GET" , requestHeaders = [("Content-type", "text/plain")] } diff --git a/src/Bio/RealWorld/ENCODE.hs b/src/Bio/RealWorld/ENCODE.hs --- a/src/Bio/RealWorld/ENCODE.hs +++ b/src/Bio/RealWorld/ENCODE.hs @@ -78,7 +78,7 @@ -- search ["chip", "sp1"] ["type=experiment"] search :: KeyWords -> IO (Either String [Value]) search kw = do - initReq <- parseUrl url + initReq <- parseRequest url let request = initReq { method = "GET" , requestHeaders = [("accept", "application/json")] } @@ -116,7 +116,7 @@ openUrl :: String -> String -> IO B.ByteString openUrl url datatype = do - initReq <- parseUrl url + initReq <- parseRequest url let request = initReq { method = "GET" , requestHeaders = [("accept", BS.pack datatype)] } diff --git a/src/Bio/RealWorld/Ensembl.hs b/src/Bio/RealWorld/Ensembl.hs --- a/src/Bio/RealWorld/Ensembl.hs +++ b/src/Bio/RealWorld/Ensembl.hs @@ -27,7 +27,7 @@ lookupHelp :: [EnsemblID] -> IO (Either String Object) lookupHelp xs = do - initReq <- parseUrl url + initReq <- parseRequest url let request = initReq { method = "POST" , requestHeaders = [("Content-type", "application/json")] , requestBody = body diff --git a/src/Bio/Seq.hs b/src/Bio/Seq.hs --- a/src/Bio/Seq.hs +++ b/src/Bio/Seq.hs @@ -1,6 +1,6 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE FlexibleContexts #-} module Bio.Seq ( -- * Alphabet @@ -19,11 +19,11 @@ , nucleotideFreq ) where -import Prelude hiding (length) import qualified Data.ByteString.Char8 as B -import qualified Data.HashMap.Strict as M -import qualified Data.HashSet as S -import Data.Char8 (toUpper) +import Data.Char8 (toUpper) +import qualified Data.HashMap.Strict as M +import qualified Data.HashSet as S +import Prelude hiding (length) -- | Alphabet defined by http://www.chem.qmul.ac.uk/iupac/ -- | Standard unambiguous alphabet @@ -103,7 +103,7 @@ 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 +-- | O(n) Reverse complementary of DNA sequence. rc :: DNA alphabet -> DNA alphabet rc (DNA s) = DNA . B.map f . B.reverse $ s where @@ -114,7 +114,7 @@ 'T' -> 'A' _ -> x --- | O(n) Compute GC content +-- | O(n) Compute GC content. gcContent :: DNA alphabet -> Double gcContent = (\(a,b) -> a / fromIntegral b) . B.foldl' f (0.0,0::Int) . toBS where @@ -133,7 +133,7 @@ _ -> x + 0.5 -- "NMKYR" in (x', n+1) --- | O(n) Compute single nucleotide frequency +-- | O(n) Compute single nucleotide frequency. nucleotideFreq :: BioSeq DNA a => DNA a -> M.HashMap Char Int nucleotideFreq dna = B.foldl' f m0 . toBS $ dna where diff --git a/src/Bio/Seq/IO.hs b/src/Bio/Seq/IO.hs --- a/src/Bio/Seq/IO.hs +++ b/src/Bio/Seq/IO.hs @@ -57,8 +57,11 @@ withGenome fl fn = bracket (openGenome fl) closeGenome fn {-# INLINE withGenome #-} -type Query = (B.ByteString, Int, Int) -- (chr, start, end), zero-based index, half-close-half-open +-- | A query is represented by a tuple: (chr, start, end) and is +-- zero-based index, half-close-half-open +type Query = (B.ByteString, Int, Int) +-- | Retrieve sequence. getSeq :: BioSeq s a => Genome -> Query -> IO (Either String (s a)) getSeq (Genome h index headerSize) (chr, start, end) = case M.lookup chr index of Just (chrStart, chrSize) -> @@ -71,6 +74,7 @@ _ -> return $ Left $ "Bio.Seq.getSeq: Cannot find " ++ show chr {-# INLINE getSeq #-} +-- | Retrieve whole chromosome. getChrom :: Genome -> B.ByteString -> IO (Either String (DNA IUPAC)) getChrom g chr = case lookup chr chrSize of Just s -> getSeq g (chr, 0, s) @@ -79,6 +83,7 @@ chrSize = getChrSizes g {-# INLINE getChrom #-} +-- | Retrieve chromosome size information. getChrSizes :: Genome -> [(B.ByteString, Int)] getChrSizes (Genome _ table _) = map (\(k, (_, l)) -> (k, l)) . M.toList $ table {-# INLINE getChrSizes #-} diff --git a/tests/Tests/Bam.hs b/tests/Tests/Bam.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/Bam.hs @@ -0,0 +1,20 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Tests.Bam (tests) where + +import Bio.Data.Bam +import Bio.Data.Bed +import Conduit +import Test.Tasty +import Test.Tasty.HUnit + +tests :: TestTree +tests = testGroup "Test: Bio.Data.Bam" + [ testCase "bamToBed" bamToBedTest + ] + +bamToBedTest :: Assertion +bamToBedTest = do + bed <- readBed' "tests/data/example.bed" + bed' <- runBam $ readBam "tests/data/example.bam" =$= bamToBed $$ sinkList + bed @=? bed' diff --git a/tests/data/example.bam b/tests/data/example.bam new file mode 100644 Binary files /dev/null and b/tests/data/example.bam differ diff --git a/tests/test.hs b/tests/test.hs --- a/tests/test.hs +++ b/tests/test.hs @@ -1,4 +1,5 @@ import qualified Tests.Bed as Bed +import qualified Tests.Bam as Bam import qualified Tests.Motif as Motif import qualified Tests.Seq as Seq import qualified Tests.GREAT as GREAT @@ -8,6 +9,7 @@ main :: IO () main = defaultMain $ testGroup "Main" [ Bed.tests + , Bam.tests , Seq.tests , Motif.tests , GREAT.tests