diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,2 +1,2 @@ -the-bioinformatician-s-toolkit -============================== +Bioinformatics Algorithms +========================= diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal --- a/bioinformatics-toolkit.cabal +++ b/bioinformatics-toolkit.cabal @@ -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 diff --git a/src/Bio/Data/Bed/Types.hs b/src/Bio/Data/Bed/Types.hs --- a/src/Bio/Data/Bed/Types.hs +++ b/src/Bio/Data/Bed/Types.hs @@ -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 diff --git a/src/Bio/Data/Fastq.hs b/src/Bio/Data/Fastq.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Data/Fastq.hs @@ -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 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 @@ -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/" diff --git a/src/Bio/RealWorld/Uniprot.hs b/src/Bio/RealWorld/Uniprot.hs new file mode 100644 --- /dev/null +++ b/src/Bio/RealWorld/Uniprot.hs @@ -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) + ] diff --git a/src/Bio/Seq.hs b/src/Bio/Seq.hs --- a/src/Bio/Seq.hs +++ b/src/Bio/Seq.hs @@ -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 diff --git a/src/Bio/Utils/Functions.hs b/src/Bio/Utils/Functions.hs --- a/src/Bio/Utils/Functions.hs +++ b/src/Bio/Utils/Functions.hs @@ -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 #-} diff --git a/tests/Tests/Tools.hs b/tests/Tests/Tools.hs --- a/tests/Tests/Tools.hs +++ b/tests/Tests/Tools.hs @@ -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] ]