diff --git a/Biobase/Types/AminoAcidSequence.hs b/Biobase/Types/AminoAcidSequence.hs new file mode 100644 --- /dev/null +++ b/Biobase/Types/AminoAcidSequence.hs @@ -0,0 +1,59 @@ + +-- | Encode the allowed amino acids in a better way. + +module Biobase.Types.AminoAcidSequence where + +import Control.DeepSeq +import Control.Lens +import Data.ByteString (ByteString) +import Data.Char (ord,chr,toUpper) +import Data.Data (Data) +import Data.Typeable (Typeable) +import GHC.Exts (IsString(..)) +import GHC.Generics (Generic) +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.UTF8 as BSU +import Test.QuickCheck (Arbitrary(..)) +import qualified Test.QuickCheck as TQ + + + +-- | A short amino acid suquence. +-- +-- It is an instance of 'Ixed' to allow @RNAseq (BS.pack "cag") ^? ix 2 == Just 'g'@. + +newtype AAseq = AAseq { _aaseq ∷ ByteString } + deriving (Data, Typeable, Generic, Eq, Ord, Read, Show) +makeLenses ''AAseq + +instance NFData AAseq + +type instance Index AAseq = Int + +type instance IxValue AAseq = Char + +instance Ixed AAseq where + ix k = aaseq . ix k . iso (chr . fromIntegral) (fromIntegral . ord) + {-# Inline ix #-} + +deriving instance Reversing AAseq + +mkAAseq ∷ ByteString → AAseq +mkAAseq = AAseq . BS.map go . BS.map toUpper + where go x | x `elem` aas = x + | otherwise = 'X' + aas ∷ String + aas = "ARNDCEQGHILKMFPSTWYVUO" + +instance IsString AAseq where + fromString = mkAAseq . BS.pack + +instance Arbitrary AAseq where + arbitrary = do + k ← TQ.choose (0,100) + xs ← TQ.vectorOf k $ TQ.elements "ARNDCEQGHILKMFPSTWYVUO" + return . AAseq $ BS.pack xs + shrink = view (to shrink) + + + diff --git a/Biobase/Types/Bitscore.hs b/Biobase/Types/Bitscore.hs --- a/Biobase/Types/Bitscore.hs +++ b/Biobase/Types/Bitscore.hs @@ -26,7 +26,7 @@ import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU -import Biobase.Types.NumericalExtremes +import Numeric.Limits @@ -48,7 +48,7 @@ instance ToJSON Bitscore instance NFData Bitscore -deriving instance NumericalExtremes Bitscore +deriving instance NumericLimits Bitscore derivingUnbox "Bitscore" [t| Bitscore -> Double |] [| getBitscore |] [| Bitscore |] @@ -58,7 +58,7 @@ -- TODO Check out the different "defaults" Infernal uses instance Default Bitscore where - def = Bitscore minLarge + def = Bitscore minFinite / 100 {-# Inline def #-} -- | Given a null model and a probability, calculate the corresponding @@ -68,7 +68,7 @@ prob2Score :: Double -> Double -> Bitscore prob2Score null x - | x==0 = minLarge + | x==0 = minFinite / 100 | otherwise = Bitscore $ log (x/null) / log 2 {-# Inline prob2Score #-} @@ -76,7 +76,7 @@ score2Prob :: Double -> Bitscore -> Double score2Prob null (Bitscore x) - | x <= minLarge = 0 + | x <= minFinite / 100 = 0 | otherwise = null * exp (x * log 2) {-# Inline score2Prob #-} diff --git a/Biobase/Types/Energy.hs b/Biobase/Types/Energy.hs --- a/Biobase/Types/Energy.hs +++ b/Biobase/Types/Energy.hs @@ -17,7 +17,8 @@ import Data.Vector.Unboxed.Deriving import GHC.Generics -import Biobase.Types.NumericalExtremes +import Numeric.Discretized +import Numeric.Limits @@ -44,35 +45,33 @@ instance ToJSON DG instance NFData DG -deriving instance NumericalExtremes DG -deriving instance NumericalEpsilon DG +deriving instance NumericLimits DG +deriving instance NumericEpsilon DG instance Default DG where - def = maxLarge + def = maxFinite / 100 {-# Inline def #-} --- | @round $ DG / 100@. - -newtype DeltaDekaGibbs = DekaG { getDekaG :: Int } - deriving (Eq,Ord,Num,Read,Show,Generic) - - - -derivingUnbox "DeltaDekaGibbs" - [t| DeltaDekaGibbs -> Int |] [| getDekaG |] [| DekaG |] - -instance Hashable DeltaDekaGibbs -instance Binary DeltaDekaGibbs -instance Serialize DeltaDekaGibbs -instance FromJSON DeltaDekaGibbs -instance ToJSON DeltaDekaGibbs -instance NFData DeltaDekaGibbs +-- | Discretized @DG@. -deriving instance NumericalExtremes DeltaDekaGibbs +newtype DDG = DDG { dDG ∷ Discretized 1 100 } + deriving (Eq,Ord,Num,Read,Show,Generic,Integral,Real,Enum) -instance Default DeltaDekaGibbs where - def = maxLarge - {-# Inline def #-} +derivingUnbox "DDG" + [t| DDG → Int |] [| getDiscretized . dDG |] [| DDG . Discretized |] +--instance Hashable DeltaDekaGibbs +--instance Binary DeltaDekaGibbs +--instance Serialize DeltaDekaGibbs +--instance FromJSON DeltaDekaGibbs +--instance ToJSON DeltaDekaGibbs +--instance NFData DeltaDekaGibbs +-- +--deriving instance NumericLimits DeltaDekaGibbs +-- +--instance Default DeltaDekaGibbs where +-- def = maxFinite `div` 100 +-- {-# Inline def #-} +-- diff --git a/Biobase/Types/Evalue.hs b/Biobase/Types/Evalue.hs --- a/Biobase/Types/Evalue.hs +++ b/Biobase/Types/Evalue.hs @@ -21,7 +21,7 @@ import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU -import Biobase.Types.NumericalExtremes +import Numeric.Limits @@ -46,17 +46,9 @@ def = Evalue 0 {-# Inline def #-} -instance NumericalExtremes Evalue where +instance NumericLimits Evalue where maxFinite = Evalue maxFinite minFinite = Evalue 0 - maxExtreme = Evalue maxExtreme - minExtreme = Evalue epsilon - maxLarge = Evalue maxLarge - minLarge = Evalue (2.2e-15) {-# Inline maxFinite #-} {-# Inline minFinite #-} - {-# Inline maxExtreme #-} - {-# Inline minExtreme #-} - {-# Inline maxLarge #-} - {-# Inline minLarge #-} diff --git a/Biobase/Types/Index.hs b/Biobase/Types/Index.hs --- a/Biobase/Types/Index.hs +++ b/Biobase/Types/Index.hs @@ -7,6 +7,8 @@ -- use it, import @Biobase.Types.Index.Type@ directly. Use @fromInt0@ to -- make clear that you count from 0 and transform to an @Index t@. I.e. -- @fromInt0 0 :: Index 1@ yields the lowest 1-base index. +-- +-- Note that internally, every lowest index starts at @0 :: Int@. module Biobase.Types.Index ( module Biobase.Types.Index @@ -20,7 +22,8 @@ import GHC.TypeLits import Text.Printf -import Biobase.Types.Index.Type +import Biobase.Types.Index.Type -- hiding (getIndex) +import qualified Biobase.Types.Index.Type as IT @@ -53,7 +56,7 @@ -- | Unsafe plus. unsafePlus :: forall t . KnownNat t => Index t -> Int -> Index t -unsafePlus i n = Index $ getIndex i + n +unsafePlus i n = Index $ IT.getIndex i + n {-# Inline unsafePlus #-} -- | Helper function that allows @subtraction@ of an 'Index' and an 'Int', @@ -66,19 +69,23 @@ -- | Delta between two 'Index' points. delta :: forall t . KnownNat t => Index t -> Index t -> Int -delta i j = abs . getIndex $ i - j +delta i j = abs . IT.getIndex $ i - j {-# Inline delta #-} -- | Unsafe minus. unsafeMinus :: forall t . KnownNat t => Index t -> Int -> Index t -unsafeMinus i n = Index $ getIndex i - n +unsafeMinus i n = Index $ IT.getIndex i - n {-# Inline unsafeMinus #-} +toInt ∷ forall t . KnownNat t ⇒ Index t → Int +{-# Inline toInt #-} +toInt i = IT.getIndex i + (fromIntegral $ natVal (Proxy ∷ Proxy t)) + -- | Return the index as an @Int@-style index that is zero-based. toInt0 :: forall t . KnownNat t => Index t -> Int -toInt0 = getIndex +toInt0 = IT.getIndex {-# Inline toInt0 #-} -- | As an index from an @Int@-style zero-based one. diff --git a/Biobase/Types/Index/Type.hs b/Biobase/Types/Index/Type.hs --- a/Biobase/Types/Index/Type.hs +++ b/Biobase/Types/Index/Type.hs @@ -8,7 +8,8 @@ import Data.Hashable (Hashable) import Data.Proxy import Data.Serialize (Serialize) -import Data.Vector.Fusion.Stream.Monadic (Step(..)) +import Data.Vector.Fusion.Stream.Monadic (Step(..), flatten) +import qualified Data.Vector.Fusion.Stream.Monadic as SM import Data.Vector.Unboxed.Deriving import GHC.Generics import GHC.TypeLits @@ -17,7 +18,6 @@ import Text.Printf import Data.PrimitiveArray.Index.Class hiding (Index) -import Data.PrimitiveArray.Vector.Compat import qualified Data.PrimitiveArray.Index.Class as PA @@ -67,19 +67,22 @@ [t| forall t . Index t -> Int |] [| getIndex |] [| Index |] instance forall t . KnownNat t => PA.Index (Index t) where - linearIndex _ _ (Index z) = z + newtype LimitType (Index t) = LtIndex Int + linearIndex (LtIndex k) (Index z) = z {-# INLINE linearIndex #-} - smallestLinearIndex (Index l) = error "still needed?" - {-# INLINE smallestLinearIndex #-} - largestLinearIndex (Index h) = h - {-# INLINE largestLinearIndex #-} - size (_) (Index h) = h + 1 + size (LtIndex h) = h + 1 {-# INLINE size #-} - inBounds (_) (Index h) (Index x) = 0<=x && x<=h + inBounds (LtIndex h) (Index x) = 0<=x && x<=h {-# INLINE inBounds #-} + zeroBound = Index 0 + {-# Inline zeroBound #-} + zeroBound' = LtIndex 0 + {-# Inline zeroBound' #-} + totalSize (LtIndex k) = [fromIntegral k] + {-# Inline totalSize #-} -instance IndexStream z => IndexStream (z:.Index t) where - streamUp (ls:.Index lf) (hs:.Index ht) = flatten mk step $ streamUp ls hs +instance (KnownNat t, IndexStream z) ⇒ IndexStream (z:.Index t) where + streamUp (ls:..LtIndex lf) (hs:..LtIndex ht) = flatten mk step $ streamUp ls hs where mk z = return (z,lf) step (z,k) | k > ht = return $ Done @@ -87,7 +90,7 @@ {-# Inline [0] mk #-} {-# Inline [0] step #-} {-# Inline streamUp #-} - streamDown (ls:.Index lf) (hs:.Index ht) = flatten mk step $ streamDown ls hs + streamDown (ls:..LtIndex lf) (hs:..LtIndex ht) = flatten mk step $ streamDown ls hs where mk z = return (z,ht) step (z,k) | k < lf = return $ Done @@ -96,7 +99,11 @@ {-# Inline [0] step #-} {-# Inline streamDown #-} -instance IndexStream (Index t) +instance (KnownNat t) ⇒ IndexStream (Index t) where + streamUp l h = SM.map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# INLINE streamUp #-} + streamDown l h = SM.map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# INLINE streamDown #-} instance Arbitrary (Index t) where arbitrary = Index <$> arbitrary diff --git a/Biobase/Types/NucleotideSequence.hs b/Biobase/Types/NucleotideSequence.hs new file mode 100644 --- /dev/null +++ b/Biobase/Types/NucleotideSequence.hs @@ -0,0 +1,217 @@ + +-- | Wrappers around biosequences. + +module Biobase.Types.NucleotideSequence where + +import Control.DeepSeq +import Control.Lens +import Data.ByteString (ByteString) +import Data.Char (ord,chr,toUpper) +import Data.Data (Data) +import Data.Typeable (Typeable) +import GHC.Exts (IsString(..)) +import GHC.Generics (Generic) +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.UTF8 as BSU +import Test.QuickCheck (Arbitrary(..)) +import qualified Test.QuickCheck as TQ + + + +-- | A sequence identifier. Just a newtype wrapped text field. Because we can +-- never know what people are up to, this is utf8-encoded. +-- +-- TODO Provide @Iso'@ for @Text@, too? +-- +-- TODO move into @Biobase.Types.SequenceID@ + +newtype SequenceID = SequenceID { _sequenceID ∷ ByteString } + deriving (Data, Typeable, Generic, Eq, Ord, Read, Show, IsString) +makeLenses ''SequenceID + +instance NFData SequenceID + +-- | Convert to a string in a unicode-aware manner. + +sequenceIDstring ∷ Iso' SequenceID String +sequenceIDstring = sequenceID . iso BSU.toString BSU.fromString +{-# Inline sequenceIDstring #-} + + + +-- | A short RNA sequence. +-- +-- It is an instance of 'Ixed' to allow @RNAseq (BS.pack "cag") ^? ix 2 == Just 'g'@. + +newtype RNAseq = RNAseq { _rnaseq ∷ ByteString } + deriving (Data, Typeable, Generic, Eq, Ord, Read, Show) +makeLenses ''RNAseq + +instance NFData RNAseq + +type instance Index RNAseq = Int + +type instance IxValue RNAseq = Char + +instance Ixed RNAseq where + ix k = rnaseq . ix k . iso (chr . fromIntegral) (fromIntegral . ord) + {-# Inline ix #-} + +deriving instance Reversing RNAseq + +mkRNAseq ∷ ByteString → RNAseq +mkRNAseq = RNAseq . BS.map go . BS.map toUpper + where go x | x `elem` acgu = x + | otherwise = 'N' + acgu ∷ String + acgu = "ACGU" + +instance IsString RNAseq where + fromString = mkRNAseq . BS.pack + +instance Arbitrary RNAseq where + arbitrary = do + k ← TQ.choose (0,100) + xs ← TQ.vectorOf k $ TQ.elements "ACGU" + return . RNAseq $ BS.pack xs + shrink = view (to shrink) + + + +-- | A short DNA sequence. +-- +-- Note everything really long should be handled by specialized libraries with +-- streaming capabilities. + +newtype DNAseq = DNAseq { _dnaseq ∷ ByteString } + deriving (Data, Typeable, Generic, Eq, Ord, Read, Show) +makeLenses ''DNAseq + +instance NFData DNAseq + +type instance Index DNAseq = Int + +type instance IxValue DNAseq = Char + +instance Ixed DNAseq where + ix k = dnaseq . ix k . iso (chr . fromIntegral) (fromIntegral . ord) + {-# Inline ix #-} + +mkDNAseq ∷ ByteString → DNAseq +mkDNAseq = DNAseq . BS.map go . BS.map toUpper + where go x | x `elem` acgt = x + | otherwise = 'N' + acgt ∷ String + acgt = "ACGT" + +instance IsString DNAseq where + fromString = mkDNAseq . BS.pack + +deriving instance Reversing DNAseq + +instance Arbitrary DNAseq where + arbitrary = do + k ← TQ.choose (0,100) + xs ← TQ.vectorOf k $ TQ.elements "ACGT" + return . DNAseq $ BS.pack xs + shrink = view (to shrink) + +-- | Simple case translation from @U@ to @T@. with upper and lower-case +-- awareness. + +rna2dna ∷ Char → Char +rna2dna = \case + 'U' → 'T' + 'u' → 't' + x → x +{-# Inline rna2dna #-} + +-- | Single character RNA complement. + +rnaComplement ∷ Char → Char +rnaComplement = \case + 'A' → 'U' + 'a' → 'u' + 'C' → 'G' + 'c' → 'g' + 'G' → 'C' + 'g' → 'c' + 'U' → 'A' + 'u' → 'a' + x → x +{-# Inline rnaComplement #-} + +-- | Simple case translation from @T@ to @U@ with upper- and lower-case +-- awareness. + +dna2rna ∷ Char → Char +dna2rna = \case + 'T' → 'U' + 't' → 'u' + x → x +{-# Inline dna2rna #-} + +-- | Single character DNA complement. + +dnaComplement ∷ Char → Char +dnaComplement = \case + 'A' → 'T' + 'a' → 't' + 'C' → 'G' + 'c' → 'g' + 'G' → 'C' + 'g' → 'c' + 'T' → 'A' + 't' → 'a' + x → x +{-# Inline dnaComplement #-} + + + +-- | Transcribes a DNA sequence into an RNA sequence. Note that 'transcribe' is +-- actually very generic. We just define its semantics to be that of +-- biomolecular transcription. +-- +-- 'transcribe' makes the assumption that, given @DNA -> RNA@, we transcribe +-- the coding strand. +-- +-- +-- @@ DNAseq "ACGT" ^. transcribe == RNAseq "ACGU" RNAseq "ACGU" ^. transcribe +-- == DNAseq "ACGT" RNAseq "ACGU" ^. from transcribe :: DNAseq == DNAseq "ACGT" +-- @@ + +class Transcribe f where + type TranscribeTo f ∷ * + transcribe ∷ Iso' f (TranscribeTo f) + +-- | Transcribe a DNA sequence into an RNA sequence. This does not @reverse@ +-- the sequence! + +instance Transcribe DNAseq where + type TranscribeTo DNAseq = RNAseq + transcribe = iso (RNAseq . BS.map dna2rna . _dnaseq) (DNAseq . BS.map rna2dna . _rnaseq) + {-# Inline transcribe #-} + +-- | Transcribe a RNA sequence into an DNA sequence. This does not @reverse@ +-- the sequence! + +instance Transcribe RNAseq where + type TranscribeTo RNAseq = DNAseq + transcribe = from transcribe + {-# Inline transcribe #-} + + + +-- | The complement of a biosequence. + +class Complement f where + complement ∷ Iso' f f + +instance Complement DNAseq where + complement = iso (DNAseq . BS.map dnaComplement . _dnaseq) (DNAseq . BS.map dnaComplement . _dnaseq) + {-# Inline complement #-} + +instance Complement RNAseq where + complement = iso (RNAseq . BS.map rnaComplement . _rnaseq) (RNAseq . BS.map rnaComplement . _rnaseq) + {-# Inline complement #-} + diff --git a/Biobase/Types/NumericalExtremes.hs b/Biobase/Types/NumericalExtremes.hs deleted file mode 100644 --- a/Biobase/Types/NumericalExtremes.hs +++ /dev/null @@ -1,76 +0,0 @@ - --- | For some values, we want to have different kind of extreme values. --- Consider a @Double@ representing an energy. We want @near infinities@ --- that do not lead to numeric problems. --- --- TODO benchmark different extremes and their interplay with algebraic --- operations. --- --- TODO consider the @ieee754@ package - -module Biobase.Types.NumericalExtremes where - - - --- | Very large and small numbers with some numerical safety to @1/0@ or --- @maxBound@ (depending on if we are @Integral@ or @RealFloat@. --- --- We have: --- --- @maxFinite >= maxExtreme >= maxLarge@ --- --- @maxLarge >= minLarge@ --- --- @minLarge >= minExtreme >= minFinite@. - -class NumericalExtremes x where - maxFinite :: x -- ^ Largest finite number - minFinite :: x -- ^ Smallest finite number - maxExtreme :: x -- ^ Around @1/ 10@ of the largest finite number - minExtreme :: x -- ^ Around @1/ 10@ of the smallest finite number - maxLarge :: x -- ^ Around @1/100@ of the largest finite number - minLarge :: x -- ^ Around @1/100@ of the smallest finite number - --- | Small numbers. - -class NumericalEpsilon x where - epsilon :: x -- ^ Smallest positive number @/= 0.0@. - - - -instance NumericalExtremes Int where - maxFinite = maxBound - minFinite = minBound - maxLarge = maxBound `div` 100 - minLarge = minBound `div` 100 - maxExtreme = maxBound `div` 10 - minExtreme = minBound `div` 10 - {-# Inline maxFinite #-} - {-# Inline minFinite #-} - {-# Inline maxExtreme #-} - {-# Inline minExtreme #-} - {-# Inline maxLarge #-} - {-# Inline minLarge #-} - - - -instance NumericalExtremes Double where - maxFinite = 1.79e+308 - minFinite = -1.79e+308 - maxExtreme = 1.79e+307 - minExtreme = -1.79e+307 - maxLarge = 1.79e+306 - minLarge = -1.79e+306 - {-# Inline maxFinite #-} - {-# Inline minFinite #-} - {-# Inline maxExtreme #-} - {-# Inline minExtreme #-} - {-# Inline maxLarge #-} - {-# Inline minLarge #-} - - - -instance NumericalEpsilon Double where - epsilon = 2.2e-16 - {-# Inline epsilon #-} - diff --git a/Biobase/Types/Odds.hs b/Biobase/Types/Odds.hs deleted file mode 100644 --- a/Biobase/Types/Odds.hs +++ /dev/null @@ -1,38 +0,0 @@ - --- | Discretized log-odds. - -module Biobase.Types.Odds where - -import Control.DeepSeq (NFData(..)) -import Data.Aeson (FromJSON,ToJSON) -import Data.Binary (Binary) -import Data.Hashable (Hashable) -import Data.Serialize (Serialize) -import Data.Vector.Unboxed.Deriving -import GHC.Generics (Generic) - - - --- | Discretized log-odds. --- --- The BLOSUM matrices, for example, store data in discretized log-odds --- form. --- --- TODO Might move up even higher into statistics modules. - -newtype DLO = DLO { getDLO :: Int } - deriving (Generic,Eq,Ord,Show,Read) - -derivingUnbox "DLO" - [t| DLO -> Int |] [| getDLO |] [| DLO |] - -instance Binary DLO -instance Serialize DLO -instance FromJSON DLO -instance ToJSON DLO -instance Hashable DLO - -instance NFData DLO where - rnf (DLO k) = rnf k - {-# Inline rnf #-} - diff --git a/Biobase/Types/ReadingFrame.hs b/Biobase/Types/ReadingFrame.hs new file mode 100644 --- /dev/null +++ b/Biobase/Types/ReadingFrame.hs @@ -0,0 +1,36 @@ + +module Biobase.Types.ReadingFrame where + +import GHC.Generics + +import Biobase.Types.Index (Index, toInt0) + + + +-- | The Reading frame. Sequence indexing starts at position 1, which starts +-- reading frame 1. Reading frame 2 and 3 start at position 2 and 3 +-- respectively. +-- +-- The reading frame should be constructed from an @Index 1@ with a smart +-- constructor to get the frame calculation right. + +newtype ReadingFrame = ReadingFrame { getReadingFrame ∷ Int } + deriving (Eq,Ord,Generic) + +nextReadingFrame ∷ ReadingFrame → ReadingFrame +{-# Inline nextReadingFrame #-} +nextReadingFrame (ReadingFrame rf) = ReadingFrame $ rf `mod` 3 + 1 + +prevReadingFrame ∷ ReadingFrame → ReadingFrame +{-# Inline prevReadingFrame #-} +prevReadingFrame (ReadingFrame rf) = ReadingFrame $ rf `mod` 3 + 2 + +-- | +-- +-- TODO should this be a type class, since we might reasonably want to +-- construct from a number of possible indices? + +fromIndex ∷ Index 1 → ReadingFrame +{-# Inline fromIndex #-} +fromIndex i = ReadingFrame $ (toInt0 i `mod` 3) + 1 + diff --git a/Biobase/Types/Sequence.hs b/Biobase/Types/Sequence.hs deleted file mode 100644 --- a/Biobase/Types/Sequence.hs +++ /dev/null @@ -1,199 +0,0 @@ - --- | Wrappers around biosequences. - -module Biobase.Types.Sequence where - -import Control.Lens -import Control.DeepSeq -import Data.ByteString (ByteString) -import Data.Char (ord,chr,toUpper) -import Data.Data (Data) -import Data.Typeable (Typeable) -import GHC.Generics (Generic) -import qualified Data.ByteString.Char8 as BS -import qualified Data.ByteString.UTF8 as BSU -import GHC.Exts (IsString(..)) - - - --- | A sequence identifier. Just a newtype wrapped text field. Because we can --- never know what people are up to, this is utf8-encoded. --- --- TODO Provide @Iso'@ for @Text@, too? - -newtype SequenceID = SequenceID { _sequenceID ∷ ByteString } - deriving (Data, Typeable, Generic, Eq, Ord, Read, Show, IsString) -makeLenses ''SequenceID - -instance NFData SequenceID - --- | Convert to a string in a unicode-aware manner. - -sequenceIDstring ∷ Iso' SequenceID String -sequenceIDstring = sequenceID . iso BSU.toString BSU.fromString -{-# Inline sequenceIDstring #-} - - - --- | A short RNA sequence. --- --- It is an instance of 'Ixed' to allow @RNAseq (BS.pack "cag") ^? ix 2 == Just 'g'@. - -newtype RNAseq = RNAseq { _rnaseq ∷ ByteString } - deriving (Data, Typeable, Generic, Eq, Ord, Read, Show) -makeLenses ''RNAseq - -instance NFData RNAseq - -type instance Index RNAseq = Int - -type instance IxValue RNAseq = Char - -instance Ixed RNAseq where - ix k = rnaseq . ix k . iso (chr . fromIntegral) (fromIntegral . ord) - {-# Inline ix #-} - -deriving instance Reversing RNAseq - -mkRNAseq ∷ ByteString → RNAseq -mkRNAseq = RNAseq . BS.map go . BS.map toUpper - where go x | x `elem` acgu = x - | otherwise = 'N' - acgu ∷ String - acgu = "ACGU" - -instance IsString RNAseq where - fromString = mkRNAseq . BS.pack - - - --- | A short DNA sequence. --- --- Note everything really long should be handled by specialized libraries with --- streaming capabilities. - -newtype DNAseq = DNAseq { _dnaseq ∷ ByteString } - deriving (Data, Typeable, Generic, Eq, Ord, Read, Show) -makeLenses ''DNAseq - -instance NFData DNAseq - -type instance Index DNAseq = Int - -type instance IxValue DNAseq = Char - -instance Ixed DNAseq where - ix k = dnaseq . ix k . iso (chr . fromIntegral) (fromIntegral . ord) - {-# Inline ix #-} - -mkDNAseq ∷ ByteString → DNAseq -mkDNAseq = DNAseq . BS.map go . BS.map toUpper - where go x | x `elem` acgt = x - | otherwise = 'N' - acgt ∷ String - acgt = "ACGT" - -instance IsString DNAseq where - fromString = mkDNAseq . BS.pack - -deriving instance Reversing DNAseq - --- | Simple case translation from @U@ to @T@. with upper and lower-case --- awareness. - -rna2dna ∷ Char → Char -rna2dna = \case - 'U' → 'T' - 'u' → 't' - x → x -{-# Inline rna2dna #-} - --- | Single character RNA complement. - -rnaComplement ∷ Char → Char -rnaComplement = \case - 'A' → 'U' - 'a' → 'u' - 'C' → 'G' - 'c' → 'g' - 'G' → 'C' - 'g' → 'c' - 'U' → 'A' - 'u' → 'a' - x → x -{-# Inline rnaComplement #-} - --- | Simple case translation from @T@ to @U@ with upper- and lower-case --- awareness. - -dna2rna ∷ Char → Char -dna2rna = \case - 'T' → 'U' - 't' → 'u' - x → x -{-# Inline dna2rna #-} - --- | Single character DNA complement. - -dnaComplement ∷ Char → Char -dnaComplement = \case - 'A' → 'T' - 'a' → 't' - 'C' → 'G' - 'c' → 'g' - 'G' → 'C' - 'g' → 'c' - 'T' → 'A' - 't' → 'a' - x → x -{-# Inline dnaComplement #-} - - - --- | Transcribes a DNA sequence into an RNA sequence. Note that 'transcribe' is --- actually very generic. We just define its semantics to be that of --- biomolecular transcription. --- --- 'transcribe' makes the assumption that, given @DNA -> RNA@, we transcribe --- the coding strand. --- --- --- @@ DNAseq "ACGT" ^. transcribe == RNAseq "ACGU" RNAseq "ACGU" ^. transcribe --- == DNAseq "ACGT" RNAseq "ACGU" ^. from transcribe :: DNAseq == DNAseq "ACGT" --- @@ - -class Transcribe f where - type TranscribeTo f ∷ * - transcribe ∷ Iso' f (TranscribeTo f) - --- | Transcribe a DNA sequence into an RNA sequence. This does not @reverse@ --- the sequence! - -instance Transcribe DNAseq where - type TranscribeTo DNAseq = RNAseq - transcribe = iso (RNAseq . BS.map dna2rna . _dnaseq) (DNAseq . BS.map rna2dna . _rnaseq) - {-# Inline transcribe #-} - --- | Transcribe a RNA sequence into an DNA sequence. This does not @reverse@ --- the sequence! - -instance Transcribe RNAseq where - type TranscribeTo RNAseq = DNAseq - transcribe = from transcribe - {-# Inline transcribe #-} - - - --- | The complement of a biosequence. - -class Complement f where - complement ∷ Iso' f f - -instance Complement DNAseq where - complement = iso (DNAseq . BS.map dnaComplement . _dnaseq) (DNAseq . BS.map dnaComplement . _dnaseq) - {-# Inline complement #-} - -instance Complement RNAseq where - complement = iso (RNAseq . BS.map rnaComplement . _rnaseq) (RNAseq . BS.map rnaComplement . _rnaseq) - {-# Inline complement #-} - diff --git a/Biobase/Types/Shape.hs b/Biobase/Types/Shape.hs new file mode 100644 --- /dev/null +++ b/Biobase/Types/Shape.hs @@ -0,0 +1,202 @@ + +-- | Shape abstractions of structures. +-- +-- Shapes do not preserve sizes of structures (say unpaired regions or stem +-- length). As such, distance measures provided here are to be used carefully! +-- +-- TODO consider how to handle the different shape levels. One option would be +-- to phantom-type everything. + +module Biobase.Types.Shape where + +import Control.DeepSeq +import Control.Lens +import Control.Monad.Error.Class +import Control.Monad (foldM,unless) +import Data.ByteString (ByteString) +import Data.Data +import Data.List (foldl1') +import Data.Monoid ((<>)) +import Data.Set (Set) +import GHC.Generics (Generic) +import qualified Data.ByteString.Char8 as BS8 +import qualified Data.List as L +import qualified Data.Set as Set + +import Data.Forest.StructuredPaired + +import qualified Biobase.Types.Structure as TS + + + +-- | Shape levels are hardcoded according to their specification. +-- +-- TODO Allow compile-time check on accepted shape levels? + +data ShapeLevel + = SL1 + | SL2 + | SL3 + | SL4 + | SL5 + deriving (Eq,Ord,Show,Read,Data,Typeable,Generic) + +instance NFData ShapeLevel + + + +-- | The type of RNA shapes. Keeps the type + +data RNAshape + = RNAshape + { _rnashapelevel ∷ !ShapeLevel + -- ^ The type of shape encoded here. + , _rnashape ∷ !ByteString + -- ^ The actual shape as a string. + } + deriving (Eq,Ord,Show,Read,Data,Typeable,Generic) +makeLenses ''RNAshape + +instance NFData RNAshape + + + +-- | Given a compactified 'SPForest', creates a shape forest of the given level. +-- +-- +-- +-- TODO needs newtyping + +shapeForest + ∷ ShapeLevel + → SPForest ByteString ByteString + → SPForest Char Char +shapeForest = preStem + where + -- | In @preStem@, we aim to close in on the next stem. @SPE@ means that we + -- reached an end in a stem. + preStem s SPE = SPE + -- | The start of a tree structure. The forest is compact, which means that + -- the element in @xs@ is, by definition, not a continuation of a stack. + preStem s (SPT _ xs _) = SPT '[' (inStem s xs) ']' + -- | + preStem s spr@(SPR rs) = inStem s spr -- = error $ "preStem/SPR " ++ show rs + -- | + preStem s (SPJ xs) + | [x] ← xs = preStem s x + -- left bulge + | [l@SPR{},x@SPT{}] ← xs = if s <= SL2 then (SPJ [SPR '_', preStem s x]) else preStem s x + -- right bulge + | [x@SPT{},r@SPR{}] ← xs = if s <= SL2 then (SPJ [preStem s x, SPR '_']) else preStem s x + | otherwise = SPJ $ map (preStem s) xs -- error $ "preStem/SPJ " ++ show xs + -- + -- | After a stem, there could be an @SPE@ element. + inStem s SPE = SPE + -- | This case happens when eradicating unstructured regions with high + -- abstraction levels. + inStem s (SPT _ xs _) = inStem s xs + inStem s (SPR rs) + | s == SL1 = SPR '_' -- = error $ "inStem / SPR " ++ show rs + | otherwise = SPE + inStem s (SPJ xs) + | [x] ← xs = error "x" + -- left bulge + | [l@SPR{},x] ← xs = if s <= SL3 then preStem s (SPJ xs) else inStem s x + -- right bulge + | [x,r@SPR{}] ← xs = if s <= SL3 then preStem s (SPJ xs) else inStem s x + -- interior loop + | [l@SPR{},x,r@SPR{}] ← xs = if s == SL5 then inStem s x else preStem s (SPJ xs) +-- | s == SL1 = error $ "inStem / SPJ " ++ show xs +-- | s == SL2 = error $ "inStem / SPJ " ++ show xs + -- multibranched loop + | otherwise = SPJ $ map (preStem s) xs + +rnass2shape lvl s = shapeForestshape lvl . shapeForest lvl . TS.compactifySPForest + . either (\e → error $ show (e,s)) id . TS.rnassSPForest $ s + +-- | turn into unit test. also reverse of the input should give reverse shape! +-- this then gives a quickcheck test, reversing the input should reverse the shape +-- +-- TODO requires generating secondary structures via @Arbitrary@. + +test lvl = shapeForestshape lvl . shapeForest lvl $ TS.compactifySPForest $ either error id $ TS.rnassSPForest $ TS.RNAss "(((((...(((..(((...))))))...(((..((.....))..)))))))).." + +{- +shapeForest SL5 = go + where + go SPE = SPE + go (SPT _ xs _) + | SPE ← xs, SPR{} ← xs, [] ← ts = SPT '[' SPE ']' + | [t] ← ts = go t + | otherwise = SPT '[' (SPJ $ map go ts) ']' + where (SPJ ys) = xs + ts = [ t | t@SPT{} ← ys ] + -- should only happen on a single unfolded structure + go (SPR _) = SPR '_' + go (SPJ xs) + | [] ← ts = SPR '_' + | [t] ← ts = go t + | otherwise = SPJ $ map go ts + where ts = [ t | t@SPT{} ← xs ] + go xs = error $ show xs ++ " should no be reached" +-} + +-- | + +shapeForestshape + ∷ ShapeLevel + → SPForest Char Char + → RNAshape +shapeForestshape lvl = RNAshape lvl . go + where + go SPE = "" + go (SPT l x r) = BS8.singleton l <> go x <> BS8.singleton r + go (SPJ xs ) = mconcat $ map go xs + go (SPR x ) = BS8.singleton x -- error "should not be reached" -- BS8.singleton x + +generateShape ∷ ShapeLevel → TS.RNAss → RNAshape +generateShape = undefined + + +-- * Distance measures on the shape string itself. + +-- | Wrapper for string-positional shapes. Intentionally chosen long name. + +data RNAshapepset = RNAshapepset { _rnashapepsetlevel ∷ ShapeLevel, _rnashapepset ∷ Set (Int,Int) } + deriving (Read,Show,Eq,Ord,Generic) +makeLenses ''RNAshapepset + +instance NFData RNAshapepset + +-- | Transform an 'RNAss' into a set of base pairs @(i,j)@. The pairs are +-- 0-based. + +rnashapePairSet + ∷ (MonadError String m) + ⇒ RNAshape + → m RNAshapepset +rnashapePairSet (RNAshape lvl s2) = do + let go (set,ks ) (i,'[') = return (set,i:ks) + go (set,i:is) (j,']') = return (Set.insert (i,j) set, is) + go (set,[] ) (j,']') = throwError $ "unequal brackets in \"" ++ BS8.unpack s2 ++ "\" at position: " ++ show j + go (set,ks ) (_,'_') = return (set,ks) + (set,ss) ← foldM go (Set.empty,[]) . L.zip [0..] $ BS8.unpack s2 + unless (null ss) $ throwError $ "unequal brackets in \"" ++ BS8.unpack s2 ++ "\" with opening bracket(s): " ++ show ss + return $ RNAshapepset lvl set +{-# Inlinable rnashapePairSet #-} + +-- | RNA pair set, but a transformation error calls @error@. + +rnassPairSet' ∷ RNAshape → RNAshapepset +rnassPairSet' = either error id . rnashapePairSet + +-- | Calculates the number of different base pairs betwwen two structures. +-- +-- TODO error out on different shape levels + +shapePairDist ∷ RNAshapepset → RNAshapepset → Int +shapePairDist (RNAshapepset lvl1 p1) (RNAshapepset lvl2 p2) = Set.size z1 + Set.size z2 + where i = Set.intersection p1 p2 + z1 = p1 `Set.difference` i + z2 = p2 `Set.difference` i + diff --git a/Biobase/Types/Strand.hs b/Biobase/Types/Strand.hs --- a/Biobase/Types/Strand.hs +++ b/Biobase/Types/Strand.hs @@ -13,51 +13,57 @@ import Data.Binary import Data.Hashable (Hashable) import Data.Serialize (Serialize) -import Data.Vector.Fusion.Stream.Monadic (Step(..)) +import Data.Vector.Fusion.Stream.Monadic (Step(..), flatten) import Data.Vector.Unboxed.Deriving import GHC.Generics import Test.QuickCheck import Text.Printf import Data.PrimitiveArray.Index.Class -import Data.PrimitiveArray.Vector.Compat +-- | Encode strand information. 'PlusStrand' is defined as the strand encoded +-- in, say, the FASTA file. 'MinusStrand' hence is the reverse complement. + newtype Strand = Strand { getStrand :: Int } deriving (Eq,Ord,Generic) instance Show Strand where - show P = "+" - show M = "-" + show PlusStrand = "PlusStrand" + show MinusStrand = "MinusStrand" instance Read Strand where readsPrec _ xs = do - ([pm],s) <- lex xs - guard $ pm `elem` ("+-PMpm" :: String) - return (go pm,s) - where go x | x `elem` ("+Pp" :: String) = P - | x `elem` ("-Mm" :: String) = M + (pm,s) <- lex xs + case pm of + "PlusStrand" → return (PlusStrand, s) + "MinusStrand" → return (MinusStrand, s) + [x] | x `elem` ("+Pp" ∷ String) → return (PlusStrand,s) + | x `elem` ("-Mm" ∷ String) → return (MinusStrand,s) + _ → [] instance Bounded Strand where - minBound = P - maxBound = M + minBound = PlusStrand + maxBound = MinusStrand instance Enum Strand where - succ P = M - succ M = error "succ M" - pred M = P - pred P = error "pred P" + succ PlusStrand = MinusStrand + succ MinusStrand = error "succ MinusStrand" + pred MinusStrand = PlusStrand + pred PlusStrand = error "pred PlusStrand" toEnum i | i>=0 && i<=1 = Strand i toEnum i = error $ "toEnum (Strand)" ++ show i fromEnum = getStrand -pattern P = Strand 0 -pattern M = Strand 1 +pattern PlusStrand = Strand 0 +pattern MinusStrand = Strand 1 -pattern Sense = P -pattern AntiSense = M +-- TODO Sense and Antisense are somewhat different +--pattern Sense = P +--pattern AntiSense = M + instance Binary Strand instance Serialize Strand instance ToJSON Strand @@ -69,19 +75,22 @@ [t| Strand -> Int |] [| getStrand |] [| Strand |] instance Index Strand where - linearIndex _ _ (Strand z) = z + newtype (LimitType Strand) = LtStrand Strand + linearIndex _ (Strand z) = z {-# INLINE linearIndex #-} - smallestLinearIndex (Strand l) = error "still needed?" - {-# INLINE smallestLinearIndex #-} - largestLinearIndex (Strand h) = h - {-# INLINE largestLinearIndex #-} - size (_) (Strand h) = h + 1 + size (LtStrand (Strand h)) = h + 1 {-# INLINE size #-} - inBounds (_) (Strand h) (Strand x) = 0<=x && x<=h + inBounds (LtStrand (Strand h)) (Strand x) = 0<=x && x<=h {-# INLINE inBounds #-} + zeroBound = Strand 0 + {-# Inline zeroBound #-} + zeroBound' = LtStrand zeroBound + {-# Inline zeroBound' #-} + totalSize (LtStrand (Strand k)) = [ fromIntegral (fromEnum k + 1) ] + {-# Inline totalSize #-} instance IndexStream z => IndexStream (z:.Strand) where - streamUp (ls:.Strand lf) (hs:.Strand ht) = flatten mk step $ streamUp ls hs + streamUp (ls:..LtStrand (Strand lf)) (hs:..LtStrand (Strand ht)) = flatten mk step $ streamUp ls hs where mk z = return (z,lf) step (z,k) | k > ht = return $ Done @@ -89,7 +98,7 @@ {-# Inline [0] mk #-} {-# Inline [0] step #-} {-# Inline streamUp #-} - streamDown (ls:.Strand lf) (hs:.Strand ht) = flatten mk step $ streamDown ls hs + streamDown (ls:..LtStrand (Strand lf)) (hs:..LtStrand (Strand ht)) = flatten mk step $ streamDown ls hs where mk z = return (z,ht) step (z,k) | k < lf = return $ Done diff --git a/Biobase/Types/Structure.hs b/Biobase/Types/Structure.hs --- a/Biobase/Types/Structure.hs +++ b/Biobase/Types/Structure.hs @@ -5,28 +5,41 @@ -- -- TODO Consider where to move each type. There are merge possibilities between -- BiobaseXNA and BiobaseTypes. +-- +-- TODO QuickCheck @Arbitrary@ for @RNAss@. module Biobase.Types.Structure where +import Control.Applicative import Control.DeepSeq import Control.Lens import Control.Monad.Error.Class import Control.Monad (foldM,unless) +import Data.Attoparsec.ByteString.Char8 +import Data.Attoparsec.Combinator +import Data.Bifunctor (second) import Data.ByteString (ByteString) import Data.Data +import Data.List (foldl1',foldl') +import Data.Monoid ((<>)) import Data.Set (Set) import GHC.Generics (Generic) import qualified Data.ByteString.Char8 as BS8 import qualified Data.List as L +import qualified Data.Set as S import qualified Data.Set as Set +import qualified Data.Vector.Unboxed as VU +import qualified Test.QuickCheck as Q +import Data.Forest.StructuredPaired + -- | Secondary structure using @()@ for paired elements, and @.@ for unpaired -- ones. It is assumed that the @()@ match up. These structures from a Monoid. newtype RNAss = RNAss { _rnass ∷ ByteString } - deriving (Eq,Ord,Show,Read,Data,Typeable,Generic,Monoid) + deriving (Eq,Ord,Show,Read,Data,Typeable,Generic,Semigroup,Monoid) makeLenses ''RNAss instance NFData RNAss @@ -100,7 +113,14 @@ verifyRNAss ss = do return ss -newtype RNApset = RNApset { _rnapset ∷ Set (Int,Int) } +-- | The set of nucleotide pairs, together with the sequence length. + +data RNApset = RNApset + { _rnapset ∷ !(Set (Int,Int)) + -- ^ the set of nucleotide pairs. + , _rnapsetSLen ∷ !Int + -- ^ length of the underlying nucleotide sequence. + } deriving (Read,Show,Eq,Ord,Generic) makeLenses ''RNApset @@ -120,19 +140,86 @@ go (set,ks ) (_,'.') = return (set,ks) (set,ss) ← foldM go (Set.empty,[]) . L.zip [0..] $ BS8.unpack s2 unless (null ss) $ throwError $ "unequal brackets in \"" ++ BS8.unpack s2 ++ "\" with opening bracket(s): " ++ show ss - return $ RNApset set + return $ RNApset set (BS8.length s2) {-# Inlinable rnassPairSet #-} +-- | Genereate a simple structured/paired forest from a secondary structure string. + +rnassSPForest + ∷ (MonadError String m) + ⇒ RNAss + → m (SPForest ByteString Char) +rnassSPForest (RNAss s2) = either throwError return $ parseOnly (manyElems <* endOfInput) s2 + where + tree = SPT <$> char '(' <*> someElems <*> char ')' "SPT" + unpaired = SPR <$> takeWhile1 (=='.') "SPR" + someElems = SPJ <$> many1 (tree <|> unpaired) "many1 SPT / SPR" + manyElems = (\case {[] → SPE; xs → SPJ xs}) <$> many (tree <|> unpaired) "many0 SPT / SPR" +{-# Inlinable rnassSPForest #-} + +-- | Compactify such an SPForest. This means that all stems are now represented +-- by a single 'SPT' data constructor. + +compactifySPForest + ∷ SPForest ByteString Char + → SPForest ByteString ByteString +compactifySPForest = go . second BS8.singleton + where go SPE = SPE + go (SPR x) = SPR x + go (SPJ xs) = SPJ (map go xs) + go (SPT l (SPJ [x]) r) = go $ SPT l x r + go (SPT l (SPT l' t r') r) = go $ SPT (l <> l') t (r' <> r) + go (SPT l t r) = SPT l (go t) r + -- | RNA pair set, but a transformation error calls @error@. rnassPairSet' ∷ RNAss → RNApset rnassPairSet' = either error id . rnassPairSet --- | Calculates the number of different base pairs betwwen two structures. +rnapsetRNAss ∷ RNApset → RNAss +rnapsetRNAss (RNApset ps l) = RNAss $ BS8.pack $ VU.toList xs + where xs = VU.replicate l '.' VU.// ls VU.// rs + ls = L.map ((,'(') . fst) $ S.toList ps + rs = L.map ((,')') . snd) $ S.toList ps +-- | Calculates the number of different base pairs between two structures. This +-- ignores the length of the underlying sequences. + pairDist ∷ RNApset → RNApset → Int -pairDist (RNApset p1) (RNApset p2) = Set.size z1 + Set.size z2 +pairDist (RNApset p1 _) (RNApset p2 _) = Set.size z1 + Set.size z2 where i = Set.intersection p1 p2 z1 = p1 `Set.difference` i z2 = p2 `Set.difference` i + + + +-- * Arbitrary instances. This only creates legal instances, but does *not* +-- take into account ViennaRNA rules like three unpaired nucleotides in the +-- hairpin. +-- +-- TODO @shrink@ is a bit more complicated, but can be done via a set of pairs. + +instance Q.Arbitrary RNApset where + arbitrary = do + -- Given left and right bounds, create pairs. + let go ∷ Int → Int → Q.Gen (Set (Int,Int)) + go l r + | l >= r = return S.empty + | otherwise = do + -- right border of stack + c ∷ Int ← Q.oneof [ Q.choose (l+1,r) -- wide jump + , Q.choose (l+1, min r $ l+20) -- short jump + ] + -- with @1..10@ stack length + z ∷ Int ← Q.choose (0,5) + let stack = S.fromList [(l+k,c-k) | k ← [0..z-1], l+k+1 < c-k] + right ← go (c+1) r + return $ S.union stack right + -- generate RNA structures between 0 and 100 nucleotides. + l ∷ Int ← Q.choose (0,199) + s ← go 0 l + return $ RNApset s (l+1) + +instance Q.Arbitrary RNAss where + arbitrary = rnapsetRNAss <$> Q.arbitrary diff --git a/Biobase/Types/Taxonomy.hs b/Biobase/Types/Taxonomy.hs --- a/Biobase/Types/Taxonomy.hs +++ b/Biobase/Types/Taxonomy.hs @@ -23,6 +23,8 @@ -- | Taxonomic classification. @Enum@ together with a final @Unknown@ is -- somewhat fishy. +-- +-- TODO What should the order be? Kingdom > Species or Kingdom < Species? data Classification = Kingdom diff --git a/BiobaseTypes.cabal b/BiobaseTypes.cabal --- a/BiobaseTypes.cabal +++ b/BiobaseTypes.cabal @@ -1,7 +1,7 @@ name: BiobaseTypes -version: 0.1.3.0 -author: Christian Hoener zu Siederdissen, 2015 - 2017 -copyright: Christian Hoener zu Siederdissen, 2015 - 2017 +version: 0.1.4.0 +author: Christian Hoener zu Siederdissen, 2015 - 2018 +copyright: Christian Hoener zu Siederdissen, 2015 - 2018 homepage: https://github.com/choener/BiobaseTypes bug-reports: https://github.com/choener/BiobaseTypes/issues maintainer: choener@bioinf.uni-leipzig.de @@ -11,7 +11,7 @@ build-type: Simple stability: experimental cabal-version: >= 1.10.0 -tested-with: GHC == 8.0.2, GHC == 8.2.1 +tested-with: GHC == 8.4.4 synopsis: Collection of types for bioinformatics description: Types used in a number of bioinformatics libraries. @@ -20,8 +20,6 @@ . * energies . - * numerics - . * biostring wrappers @@ -35,6 +33,7 @@ library build-depends: base >= 4.7 && < 5.0 , aeson >= 0.8 + , attoparsec >= 0.13 , binary >= 0.7 , bytestring , cereal >= 0.4 @@ -58,9 +57,12 @@ , vector-th-unbox >= 0.2 -- , bimaps == 0.1.0.* - , PrimitiveArray == 0.8.0.* + , ForestStructures == 0.0.1.* + , PrimitiveArray == 0.9.0.* + , SciBaseTypes == 0.0.0.* exposed-modules: Biobase.Types.Accession + Biobase.Types.AminoAcidSequence Biobase.Types.Bitscore Biobase.Types.Energy Biobase.Types.Evalue @@ -68,9 +70,9 @@ Biobase.Types.Index.Type Biobase.Types.Names Biobase.Types.Names.Internal - Biobase.Types.NumericalExtremes - Biobase.Types.Odds - Biobase.Types.Sequence + Biobase.Types.NucleotideSequence + Biobase.Types.ReadingFrame + Biobase.Types.Shape Biobase.Types.Strand Biobase.Types.Structure Biobase.Types.Taxonomy @@ -88,11 +90,14 @@ , MultiParamTypeClasses , OverloadedStrings , PatternSynonyms + , PolyKinds , ScopedTypeVariables , StandaloneDeriving , TemplateHaskell + , TypeApplications , TypeFamilies , TypeOperators + , TupleSections , UnicodeSyntax ghc-options: -O2 -funbox-strict-fields @@ -112,7 +117,10 @@ Haskell2010 default-extensions: ScopedTypeVariables , TemplateHaskell + , UnicodeSyntax build-depends: base + , bytestring + , lens , QuickCheck , tasty >= 0.11 , tasty-quickcheck >= 0.8 diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,9 @@ +0.1.4.0 +------- + +- changes to indexing and others +- some changes that probably require a bump + 0.1.3.0 ------- diff --git a/tests/properties.hs b/tests/properties.hs --- a/tests/properties.hs +++ b/tests/properties.hs @@ -1,6 +1,9 @@ module Main where +import Control.Lens +import Debug.Trace +import qualified Data.ByteString.Char8 as BS8 import Test.QuickCheck.Modifiers import Test.QuickCheck.Property import Test.Tasty @@ -8,7 +11,9 @@ import Test.Tasty.TH import Biobase.Types.Bitscore -import Biobase.Types.NumericalExtremes +import Biobase.Types.NucleotideSequence +import Biobase.Types.Shape +import Biobase.Types.Structure @@ -18,6 +23,49 @@ --prop_ScoreProb (Positive null) x = Bitscore x ~= prob2Score null (score2Prob null $ Bitscore x) + + +-- * sequence properties + +-- complement twice + +prop_complement_twice_DNA (dna ∷ DNAseq) = dna == dna^.complement.complement + +prop_complement_twice_RNA (rna ∷ RNAseq) = rna == rna^.complement.complement + +prop_transcribe_twice_DNA (dna ∷ DNAseq) = dna == dna^.transcribe.transcribe + +--prop_transcribe_twice_DNA (rna ∷ RNAseq) = rna == rna^.transcribe.transcribe + +-- * shape properties + +-- ** unit tests for known rna secondary structures + +-- ** quickcheck + +-- | reversing a secondary structure means reversing the shape + +prop_StructureShape_5_Reverse = fun_StructureShape_k_Reverse SL5 +prop_StructureShape_4_Reverse = fun_StructureShape_k_Reverse SL4 +prop_StructureShape_3_Reverse = fun_StructureShape_k_Reverse SL3 +prop_StructureShape_2_Reverse = fun_StructureShape_k_Reverse SL2 +prop_StructureShape_1_Reverse = fun_StructureShape_k_Reverse SL1 + +fun_StructureShape_k_Reverse lvl rnass@(RNAss s2) + | shp == fshp = True + | otherwise = traceShow (s2,shp,rshp,fshp) False + where shp = rnass2shape lvl rnass + rshp = rnass2shape lvl $ RNAss $ BS8.map flp $ BS8.reverse s2 + fshp = over rnashape (BS8.map flp . BS8.reverse) rshp + flp '(' = ')' + flp ')' = '(' + flp '[' = ']' + flp ']' = '[' + flp x = x + + + +-- * generic stuff a ~= b = abs (b-a) <= 10e-6