BiobaseTypes 0.1.3.0 → 0.2.1.0
raw patch · 24 files changed
Files
- Biobase/Types/BioSequence.hs +409/−0
- Biobase/Types/Bitscore.hs +18/−6
- Biobase/Types/Codon.hs +20/−0
- Biobase/Types/Energy.hs +38/−21
- Biobase/Types/Evalue.hs +4/−10
- Biobase/Types/Index.hs +33/−19
- Biobase/Types/Index/Type.hs +30/−17
- Biobase/Types/Location.hs +253/−0
- Biobase/Types/NumericalExtremes.hs +0/−76
- Biobase/Types/Odds.hs +0/−38
- Biobase/Types/Position.hs +163/−0
- Biobase/Types/ReadingFrame.hs +54/−0
- Biobase/Types/Sequence.hs +0/−199
- Biobase/Types/Shape.hs +202/−0
- Biobase/Types/Strand.hs +72/−35
- Biobase/Types/Structure.hs +92/−5
- Biobase/Types/Taxonomy.hs +2/−0
- BiobaseTypes.cabal +54/−32
- DP/Backtraced/BioSequence.hs +3/−0
- DP/Backtraced/Codon.hs +46/−0
- LICENSE +1/−1
- README.md +2/−1
- changelog.md +22/−0
- tests/properties.hs +111/−3
+ Biobase/Types/BioSequence.hs view
@@ -0,0 +1,409 @@++-- | Abstraction over bio sequences encoded as one-ascii character as one+-- symbol. We phantom-type the exact bio-sequence type and provide type classes+-- that act on known types.+--+-- Unknown bio sequences should be tagged with @Void@.+--+-- TODO give (lens) usage examples++module Biobase.Types.BioSequence where++import Control.DeepSeq+import Control.Lens+import Data.ByteString.Char8 (ByteString)+import Data.Char (ord,chr,toUpper)+import Data.Data (Data)+import Data.Hashable+import Data.Typeable (Typeable)+import Data.Void+import GHC.Exts (IsString(..))+import GHC.Generics (Generic)+import qualified Data.ByteString.Char8 as BS+import qualified Data.ByteString.UTF8 as BSU+import qualified Streaming.Prelude as SP+import qualified Streaming as S+import qualified Streaming.Internal as SI+import qualified Test.QuickCheck as TQ+import Test.QuickCheck (Arbitrary(..))+import Data.Coerce+import Debug.Trace++import Biobase.Types.Strand+import qualified Biobase.Types.Index as BTI+import Data.Info++++-- * Lens operations on biosequences++{-+class BioSeqLenses b where+ -- | Lens into the first @k@ characters.+ bsTake :: Int -> Lens' b b+ -- | Lens into the last @k@ characters+ bsTakeEnd :: Int -> Lens' b b+ -- | Lens into all but the first @k@ characters+ bsDrop :: Int -> Lens' b b+ -- | Lens into all but the last @k@ characters+ bsDropEnd :: Int -> Lens' b b+ -- | Lens that splits at a position+ bsSplitAt :: Int -> Lens' b (b,b)+ -- | length of this biosequence+ bsLength :: Getter b Int+-}++-- * Sequence identifiers++newtype SequenceIdentifier (which :: k) = SequenceIdentifier { _sequenceIdentifier :: ByteString }+ deriving stock (Data, Typeable, Generic, Eq, Ord, Read, Show)+makeWrapped ''SequenceIdentifier+makePrisms ''SequenceIdentifier++instance NFData (SequenceIdentifier w)++instance IsString (SequenceIdentifier w) where+ fromString = SequenceIdentifier . BSU.fromString++++-- * Bio-Sequences++data RNA++data DNA++data XNA++data AA++++-- |+-- TODO provide extended annotation information on biosequences, too!++newtype BioSequence (which :: k) = BioSequence {_bioSequence :: ByteString}+ deriving stock (Data, Typeable, Generic, Eq, Ord, Read, Show)+ deriving newtype (Semigroup)+makeWrapped ''BioSequence+makePrisms ''BioSequence+makeLenses ''BioSequence++instance Hashable (BioSequence (which :: k))++instance NFData (BioSequence w)++type instance Index (BioSequence w) = Int++type instance IxValue (BioSequence w) = Char++instance Ixed (BioSequence w) where+ ix k = _BioSequence . ix k . iso (chr . fromIntegral) (fromIntegral . ord)+ {-# Inline ix #-}++deriving newtype instance Reversing (BioSequence w)++instance IsString (BioSequence Void) where+ fromString = BioSequence . BS.pack++instance Info (BioSequence w) where+ info (BioSequence s)+ | BS.length s <= 18 = BS.unpack s+ | otherwise = BS.unpack h ++ ".." ++ BS.unpack l+ where (h,tl) = BS.splitAt 9 s+ (_,l ) = BS.splitAt (BS.length tl-9) tl++{-+instance BioSeqLenses (BioSequence w) where+ {-# Inline bsTake #-}+ bsTake k = lens (over _BioSequence (BS.take k)) (\old new -> new <> over _BioSequence (BS.drop k) old)+ {-# Inline bsTakeEnd #-}+ bsTakeEnd k = lens (over _BioSequence (\s -> BS.drop (BS.length s -k) s)) (\old new -> over _BioSequence (\s -> BS.take (BS.length s-k) s) old <> new)+ {-# Inline bsLength #-}+ bsLength = _BioSequence.to BS.length+ {-# Inline bsDrop #-}+ bsDrop k = lens (over _BioSequence (BS.drop k)) (\old new -> over _BioSequence (BS.take k) old <> new)+ {-# Inline bsDropEnd #-}+ bsDropEnd k = lens (over _BioSequence (\s -> BS.take (BS.length s -k) s)) (\old new -> over _BioSequence (\s -> BS.take (BS.length s-k) s) old <> new)+ {-# Inline bsSplitAt #-}+ bsSplitAt k = lens (\b -> (view (bsTake k) b, view (bsDrop k) b)) (\old (h,t) -> h <> t)+-}++++-- * RNA++-- |+--+-- TODO write that converts explicitly++mkRNAseq :: ByteString -> BioSequence RNA+mkRNAseq = BioSequence . BS.map go . BS.map toUpper+ where go x | x `elem` acgu = x+ | otherwise = 'N'+ acgu :: String+ acgu = "ACGU"++instance IsString (BioSequence RNA) where+ fromString = mkRNAseq . BS.pack++instance Arbitrary (BioSequence RNA) where+ arbitrary = do+ k ← TQ.choose (0,30)+ xs ← TQ.vectorOf k $ TQ.elements "ACGU"+ return . BioSequence $ BS.pack xs+ shrink = shrinkBioSequence++shrinkBioSequence (BioSequence b) = fmap BioSequence+ [ let (l,BS.drop 1 -> r) = BS.splitAt k b+ in BS.append l r | k <- [0 .. BS.length b -1] ]+++-- * DNA++mkDNAseq :: ByteString -> (BioSequence DNA)+mkDNAseq = BioSequence . BS.map go . BS.map toUpper+ where go x | x `elem` acgt = x+ | otherwise = 'N'+ acgt :: String+ acgt = "ACGT"++instance IsString (BioSequence DNA) where+ fromString = mkDNAseq . BS.pack++instance Arbitrary (BioSequence DNA) where+ arbitrary = do+ k ← TQ.choose (0,100)+ xs ← TQ.vectorOf k $ TQ.elements "ACGT"+ return . BioSequence $ BS.pack xs+ shrink = view (to shrink)++++-- * XNA++mkXNAseq :: ByteString -> (BioSequence XNA)+mkXNAseq = BioSequence . BS.map go . BS.map toUpper+ where go x | x `elem` acgtu = x+ | otherwise = 'N'+ acgtu :: String+ acgtu = "ACGTU"++instance IsString (BioSequence XNA) where+ fromString = mkXNAseq . BS.pack++instance Arbitrary (BioSequence XNA) where+ arbitrary = do+ k ← TQ.choose (0,100)+ xs ← TQ.vectorOf k $ TQ.elements "ACGTU"+ return . BioSequence $ BS.pack xs+ shrink = view (to shrink)++++-- * Amino acid sequences++mkAAseq :: ByteString -> (BioSequence AA)+mkAAseq = BioSequence . BS.map go . BS.map toUpper+ where go x | x `elem` aas = x+ | otherwise = 'X'+ aas :: String+ aas = "ARNDCEQGHILKMFPSTWYVUO"++instance IsString (BioSequence AA) where+ fromString = mkAAseq . BS.pack++instance Arbitrary (BioSequence AA) where+ arbitrary = do+ k ← TQ.choose (0,100)+ xs ← TQ.vectorOf k $ TQ.elements "ARNDCEQGHILKMFPSTWYVUO"+ return . BioSequence $ BS.pack xs+ shrink = view (to shrink)++++{-++-- * A window into a longer sequence with prefix/suffix information.++-- | Phantom-typed over two types, the type @w@ of the identifier, which can be+-- descriptive ("FirstInput") and the second type, identifying what kind of+-- sequence types we are dealing with. Finally, the third type provides+-- location information and should be location or streamed location.++data BioSequenceWindow w ty loc = BioSequenceWindow+ { _bswIdentifier :: !(SequenceIdentifier w)+ -- ^ Identifier for this window. Typically some fasta identifier+ , _bswPrefix :: !(BioSequence ty)+ , _bswInfix :: !(BioSequence ty)+ , _bswSuffix :: !(BioSequence ty)+ , _bswInfixLocation :: !loc+ -- ^ Location of the infix sequence+ }+ deriving (Data, Typeable, Generic, Eq, Ord, Read, Show)+makeLenses ''BioSequenceWindow++-- | Lens into the full sequence. May not change the sequence length++bswSequence :: Lens (BioSequenceWindow w ty loc) (BioSequenceWindow w ty' loc) (BioSequence ty) (BioSequence ty')+{-# Inlinable bswSequence #-}+bswSequence = lens (\w -> _bswPrefix w <> _bswInfix w <> _bswSuffix w)+ (\w bs -> let (p,is) = bs^.bsSplitAt (w^.bswPrefix.bsLength)+ (i,s ) = is^.bsSplitAt (w^.bswInfix.bsLength)+ in w { _bswPrefix = p, _bswInfix = i, _bswSuffix = s } )++-- | Get the position of the whole sequence++bswLocation :: ModifyLocation loc => Getter (BioSequenceWindow w ty loc) loc+{-# Inlinable bswLocation #-}+bswLocation = to $ \w -> locMoveLeftEnd (w^.bswPrefix.bsLength.to negate)+ . locMoveRightEnd (w^.bswSuffix.bsLength) $ w^.bswInfixLocation++bswRetagW :: BioSequenceWindow w ty loc -> BioSequenceWindow v ty loc+{-# Inlinable bswRetagW #-}+bswRetagW = over bswIdentifier coerce++instance NFData loc => NFData (BioSequenceWindow w ty loc)++instance (Reversing loc) => Reversing (BioSequenceWindow w ty loc) where+ {-# Inlinable reversing #-}+ reversing bsw = bsw+ & bswPrefix .~ (bsw^.bswSuffix.reversed)+ & bswSuffix .~ (bsw^.bswPrefix.reversed)+ & bswInfix .~ (bsw^.bswInfix.reversed)+ & bswInfixLocation .~ (bsw^.bswInfixLocation.reversed)++++-- | Provides an informative string indicating the current window being worked on. Requires length+-- of pretty string requested. Not for computers, but for logging what is being worked on. Should be+-- one line at most, not produce line breaks.+--+-- @...PFX [Start] IFX...IFX [End] SFX ...@+--+-- TODO possibly be better as a @Doc@ for prettier printing.++instance Info (BioSequenceWindow w ty loc) where+ info bsw = "todo: info bsw"++-}++++-- * DNA/RNA++-- | 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.+-- <http://hyperphysics.phy-astr.gsu.edu/hbase/Organic/transcription.html>+--+-- @@ 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 (BioSequence DNA) where+ type TranscribeTo (BioSequence DNA) = (BioSequence RNA)+ transcribe = iso (over _BioSequence (BS.map dna2rna)) (over _BioSequence (BS.map rna2dna))+ {-# Inline transcribe #-}++-- | Transcribe a RNA sequence into an DNA sequence. This does not @reverse@+-- the sequence!++instance Transcribe (BioSequence RNA) where+ type TranscribeTo (BioSequence RNA) = (BioSequence DNA)+ transcribe = from transcribe+ {-# Inline transcribe #-}++++-- | The complement of a biosequence.++class Complement f where+ complement :: Iso' f f++instance Complement (BioSequence DNA) where+ {-# Inline complement #-}+ complement = let f = (over _BioSequence (BS.map dnaComplement))+ {-# Inline f #-}+ in iso f f++instance Complement (BioSequence RNA) where+ {-# Inline complement #-}+ complement = let f = (over _BioSequence (BS.map rnaComplement))+ {-# Inline f #-}+ in iso f f++{-+instance (Complement (BioSequence ty)) => Complement (BioSequenceWindow w ty k) where+ {-# Inline complement #-}+ complement = let f = over bswPrefix (view complement) . over bswInfix (view complement) . over bswSuffix (view complement)+ {-# Inline f #-}+ in iso f f+-}++reverseComplement :: (Complement f, Reversing f) => Iso' f f+{-# Inline reverseComplement #-}+reverseComplement = reversed . complement+
Biobase/Types/Bitscore.hs view
@@ -26,7 +26,8 @@ import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU -import Biobase.Types.NumericalExtremes+import Algebra.Structure.Semiring+import Numeric.Limits @@ -39,8 +40,19 @@ -- Infernal users guide, p.42: log-odds score in log_2 (aka bits). newtype Bitscore = Bitscore { getBitscore :: Double }- deriving (Eq,Ord,Read,Show,Num,Fractional,Generic)+ deriving stock (Eq,Ord,Read,Show,Generic)+ deriving newtype (Num,Fractional) +instance Semiring Bitscore where+ plus = (+)+ times = (*)+ zero = 0+ one = 1+ {-# Inline plus #-}+ {-# Inline times #-}+ {-# Inline zero #-}+ {-# Inline one #-}+ instance Binary Bitscore instance FromJSON Bitscore instance Hashable Bitscore@@ -48,7 +60,7 @@ instance ToJSON Bitscore instance NFData Bitscore -deriving instance NumericalExtremes Bitscore+deriving newtype instance NumericLimits Bitscore derivingUnbox "Bitscore" [t| Bitscore -> Double |] [| getBitscore |] [| Bitscore |]@@ -58,7 +70,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 +80,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 +88,7 @@ score2Prob :: Double -> Bitscore -> Double score2Prob null (Bitscore x)- | x <= minLarge = 0+ | x <= minFinite / 100 = 0 | otherwise = null * exp (x * log 2) {-# Inline score2Prob #-}
+ Biobase/Types/Codon.hs view
@@ -0,0 +1,20 @@++module Biobase.Types.Codon where++import Control.Lens+import GHC.Generics (Generic)++++-- | A single codon.+--+-- TODO needs to go into its own place++data Codon c = Codon !c !c !c+ deriving (Eq,Ord,Read,Show,Generic,Functor,Foldable,Traversable)++instance Field1 (Codon c) (Codon c) c c+instance Field2 (Codon c) (Codon c) c c+instance Field3 (Codon c) (Codon c) c c+instance Each (Codon c) (Codon c') c c'+
Biobase/Types/Energy.hs view
@@ -13,11 +13,14 @@ import Data.Data import Data.Default import Data.Hashable+import GHC.Real import Data.Serialize (Serialize) import Data.Vector.Unboxed.Deriving import GHC.Generics -import Biobase.Types.NumericalExtremes+import Algebra.Structure.Semiring+import Numeric.Discretized+import Numeric.Limits @@ -44,35 +47,49 @@ 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)-+-- | Discretized @DG@. +newtype DDG = DDG { dDG ∷ Discretized (1 :% 100) }+ deriving (Eq,Ord,Num,Read,Generic,Real,Enum) -derivingUnbox "DeltaDekaGibbs"- [t| DeltaDekaGibbs -> Int |] [| getDekaG |] [| DekaG |]+instance Show DDG where+ show (DDG e) = show e -instance Hashable DeltaDekaGibbs-instance Binary DeltaDekaGibbs-instance Serialize DeltaDekaGibbs-instance FromJSON DeltaDekaGibbs-instance ToJSON DeltaDekaGibbs-instance NFData DeltaDekaGibbs+ddg2Int :: DDG -> Int+ddg2Int (DDG (Discretized e)) = e -deriving instance NumericalExtremes DeltaDekaGibbs+derivingUnbox "DDG"+ [t| DDG -> Int |] [| getDiscretized . dDG |] [| DDG . Discretized |] -instance Default DeltaDekaGibbs where- def = maxLarge- {-# Inline def #-}+instance Semiring DDG where+ plus (DDG x) (DDG y) = DDG $ min x y+ times (DDG x) (DDG y) = DDG $ x `plus` y+ zero = DDG maxFinite+ one = DDG zero+ {-# Inline plus #-}+ {-# Inline times #-}+ {-# Inline zero #-}+ {-# Inline one #-} +--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 #-}+--
Biobase/Types/Evalue.hs view
@@ -8,6 +8,7 @@ module Biobase.Types.Evalue where import Control.DeepSeq+import Control.Lens import Data.Aeson import Data.Binary import Data.Default@@ -21,7 +22,7 @@ import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU -import Biobase.Types.NumericalExtremes+import Numeric.Limits @@ -29,6 +30,7 @@ newtype Evalue = Evalue { getEvalue :: Double } deriving (Eq,Ord,Read,Show,Num,Generic)+makeWrapped ''Evalue instance Binary Evalue instance FromJSON Evalue@@ -46,17 +48,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 #-}
Biobase/Types/Index.hs view
@@ -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@@ -16,11 +18,13 @@ , Index ) where +import Data.Coerce import Data.Proxy 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 @@ -37,11 +41,13 @@ -- as @0 :: Index 0@ gives @1 :: Index 1@ for example. I.e. valid indices -- become valid indices. -reIndex :: forall n m . (KnownNat n, KnownNat m) => Index n -> Index m-reIndex (Index i) = Index $ i - n + m- where n = fromIntegral $ natVal (Proxy :: Proxy n)- m = fromIntegral $ natVal (Proxy :: Proxy m)+reIndex ∷ Index n → Index m {-# Inline reIndex #-}+reIndex = coerce+--reIndex :: forall n m . (KnownNat n, KnownNat m) => Index n -> Index m+--reIndex (Index i) = Index $ i - n + m+-- where n = fromIntegral $ natVal (Proxy :: Proxy n)+-- m = fromIntegral $ natVal (Proxy :: Proxy m) -- | Helper function that allows @addition@ of an 'Index' and an 'Int', -- with the 'Int' on the right.@@ -50,36 +56,44 @@ (+.) i n = checkIndex $ unsafePlus i n {-# Inline (+.) #-} --- | Unsafe plus.--unsafePlus :: forall t . KnownNat t => Index t -> Int -> Index t-unsafePlus i n = Index $ getIndex i + n-{-# Inline unsafePlus #-}- -- | Helper function that allows @subtraction@ of an 'Index' and an 'Int', -- with the 'Int' on the right. (-.) :: forall t . KnownNat t => Index t -> Int -> Index t-(-.) i n = checkIndex $ unsafeMinus i n+(-.) i n = checkIndex $ unsafePlus i (negate n) {-# Inline (-.) #-} +-- | Unsafe plus.++unsafePlus :: forall t . KnownNat t => Index t -> Int -> Index t+unsafePlus i n = Index $ IT.getIndex i + n+{-# Inline unsafePlus #-}+ -- | Delta between two 'Index' points. delta :: forall t . KnownNat t => Index t -> Index t -> Int-delta i j = abs . getIndex $ i - j+delta (Index i) (Index j) = abs $ i - j {-# Inline delta #-} --- | Unsafe minus.--unsafeMinus :: forall t . KnownNat t => Index t -> Int -> Index t-unsafeMinus i n = Index $ 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 #-}++-- | Return the index as an @Int@-style index that is one-based.++toInt1 ∷ forall t . KnownNat t ⇒ Index t → Int+{-# Inline toInt1 #-}+toInt1 = (+1) . toInt0++fromInt1 ∷ forall t . KnownNat t ⇒ Int → Index t+{-# Inline fromInt1 #-}+fromInt1 = fromInt0 . (subtract 1) -- | As an index from an @Int@-style zero-based one. --
Biobase/Types/Index/Type.hs view
@@ -5,19 +5,21 @@ import Control.DeepSeq import Data.Aeson import Data.Binary+import Data.Data (Data) import Data.Hashable (Hashable) import Data.Proxy import Data.Serialize (Serialize)-import Data.Vector.Fusion.Stream.Monadic (Step(..))+import Data.Typeable (Typeable)+import Data.Vector.Fusion.Stream.Monadic (Step(..), flatten) import Data.Vector.Unboxed.Deriving import GHC.Generics import GHC.TypeLits import qualified Data.Ix as Ix+import qualified Data.Vector.Fusion.Stream.Monadic as SM import Test.QuickCheck import Text.Printf import Data.PrimitiveArray.Index.Class hiding (Index)-import Data.PrimitiveArray.Vector.Compat import qualified Data.PrimitiveArray.Index.Class as PA @@ -25,7 +27,7 @@ -- | A linear @Int@-based index type. newtype Index (t :: Nat) = Index { getIndex :: Int }- deriving (Show,Read,Eq,Ord,Generic,Ix.Ix)+ deriving (Show,Read,Eq,Ord,Generic,Ix.Ix,Data,Typeable) -- | Turn an 'Int' into an 'Index' safely. @@ -44,9 +46,9 @@ {-# Inline maybeIndex #-} instance KnownNat t => Num (Index t) where- Index a + Index b = error "not implemented, use (+.)" -- index $ a + b- Index a - Index b = error "not implemented, use (-.)" -- index $ a - b- Index a * Index b = error "not implemented" -- index $ a * b+ Index a + Index b = error $ show (" Index.(+) not implemented, use (+.)",a,b) -- index $ a + b+ Index a - Index b = error $ show (" Index.(-) not implemented, use (-.)",a,b) -- index $ a - b+ Index a * Index b = error $ show (" Index.(*) not implemented", a,b) -- index $ a * b negate = error "Indices are natural numbers" abs = id signum = index . signum . getIndex@@ -67,19 +69,26 @@ [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 #-}+ fromLinearIndex _ = Index+ {-# Inline [0] fromLinearIndex #-}+ showBound (LtIndex k) = ["LtIndex " ++ show k]+ showIndex (Index k) = ["Index " ++ show k] -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 +96,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 +105,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
+ Biobase/Types/Location.hs view
@@ -0,0 +1,253 @@++-- | Annotate the genomic @Location@ of features or elements. A @Location@ is+-- always contiguous, using strand, 0-based position, and length.+-- Transformation to different systems of annotation is made possible.++module Biobase.Types.Location where++import Control.DeepSeq+import Control.Lens hiding (Index, index)+import Data.Coerce+import Data.Data+import Data.Data.Lens+import GHC.Generics (Generic)+import GHC.TypeNats+import Prelude hiding (length)+import qualified Data.ByteString as BS+import qualified Streaming.Internal as SI+import qualified Streaming.Prelude as SP+import Text.Printf++import Biobase.Types.BioSequence+import Biobase.Types.Index+import Biobase.Types.Position+import Biobase.Types.Strand+import Data.Info+++++-- | Operations on locations.++class ModifyLocation posTy seqTy where+ -- | Append to the left.+ locAppendLeft :: seqTy -> Location i posTy seqTy -> Location i posTy seqTy+ -- | Append to the right.+ locAppendRight :: seqTy -> Location i posTy seqTy -> Location i posTy seqTy+ -- | Split a location.+ locSplitAt :: Int -> Location i posTy seqTy -> (Location i posTy seqTy, Location i posTy seqTy)+ -- | Length of location+ locLength :: Location i posTy seqTy -> Int++locTake k = fst . locSplitAt k++locTakeEnd k loc = let l = locLength loc in snd $ locSplitAt (l-k) loc++locDrop k = snd . locSplitAt k++locDropEnd k loc = let l = locLength loc in fst $ locSplitAt (l-k) loc++locSplitEndAt k loc = let l = locLength loc in locSplitAt (l-k) loc++++data Location ident posTy seqTy = Location+ { _locIdentifier :: !(SequenceIdentifier ident)+ , _locPosition :: !posTy+ , _locSequence :: !seqTy+ }+ deriving stock (Show,Data,Typeable,Generic)+makeLenses ''Location++instance (NFData p, NFData s) => NFData (Location i p s)++retagLocation :: Location i posTy seqTy -> Location j posTy seqTy+{-# Inline retagLocation #-}+retagLocation = over locIdentifier coerce++instance ModifyLocation FwdPosition (BioSequence w) where+ {-# Inline locAppendLeft #-}+ locAppendLeft s loc = let l = s^._BioSequence.to BS.length+ in loc & locSequence %~ (s <>) & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (-. l) else p)+ {-# Inline locAppendRight #-}+ locAppendRight s loc = let l = s^._BioSequence.to BS.length+ in loc & locSequence %~ (<> s) & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~ (-. l) else p)+ {-# Inline locSplitAt #-}+ locSplitAt k loc =+ let (h',t') = loc^.locSequence._BioSequence.to (BS.splitAt k)+ hl = BS.length h' ; tl = BS.length t'+ h = loc & locSequence._BioSequence .~ h' & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~ (+. tl) else p)+ t = loc & locSequence._BioSequence .~ t' & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (+. hl) else p)+ in (h,t)+ {-# Inline locLength #-}+ locLength = view (locSequence._BioSequence.to BS.length)++instance ModifyLocation FwdPosition Int where+ {-# Inline locAppendLeft #-}+ locAppendLeft k' loc = let k = max 0 $ min (loc^.locPosition.fwdStart.to toInt0) k' in+ loc & locSequence %~ (+ k) & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (-. k) else p)+ {-# Inline locAppendRight #-}+ locAppendRight k' loc = let k = max 0 $ min (loc^.locPosition.fwdStart.to toInt0) k' in+ loc & locSequence %~ (+ k) & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~ (-. k) else p)+ {-# Inline locSplitAt #-}+ locSplitAt k loc =+ let h' = max 0 . min k $ locLength loc+ t' = locLength loc - h'+ h = loc & locSequence .~ h' & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~(+. t') else p)+ t = loc & locSequence .~ t' & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (+. h') else p)+ in (h,t)+ {-# Inline locLength #-}+ locLength = view locSequence++instance Reversing (Location i FwdPosition (BioSequence w)) where+ {-# Inline reversing #-}+ reversing = over (locSequence._BioSequence) BS.reverse . over (locPosition) reversing++instance Complement (BioSequence w) => Complement (Location i FwdPosition (BioSequence w)) where+ {-# Inline complement #-}+ complement = iso f f+ where f = over locSequence (view complement)++instance (Info (BioSequence w)) => Info (Location i FwdPosition (BioSequence w)) where+ info loc = printf "%s %s %s" (loc^.locIdentifier^.to show) (show $ loc^.locPosition) (loc^.locSequence.to info)++-- | Will extract a substring for a given biosequence. It is allowed to hand in partially or not at+-- all overlapping locational information. This will yield empty resulting locations.+--+-- This will convert the @FwdPosition@ strand, which in turn allows dealing with reverse-complement+-- searches.+--+-- @+-- 0123456789+-- 3.3+-- @++subLocation :: Location i FwdPosition (BioSequence w) -> (FwdPosition, Int) -> Location i FwdPosition (BioSequence w)+{-# Inline subLocation #-}+subLocation s (p',l)+ | ss==PlusStrand = locTake l $ locDrop d s+ | ss==MinusStrand = locTakeEnd l $ locDropEnd d s+ where ss = s^.locPosition.fwdStrand+ p = if ss == p'^.fwdStrand then p' else reversing p'+ d = delta (s^.locPosition.fwdStart) (p^.fwdStart)++data PIS i p s = PIS+ { _pisPrefix :: Maybe (Location i p s)+ , _pisInfix :: !(Location i p s)+ , _pisSuffix :: Maybe (Location i p s)+ }+ deriving stock (Show, Data)+makeLenses ''PIS++pis ifx = PIS Nothing ifx Nothing++retagPis :: PIS i p s -> PIS j p s+retagPis (PIS p i s) = PIS (fmap retagLocation p) (retagLocation i) (fmap retagLocation s)++-- | Given a @PIS@, this will return the @substring@ indicated by the location in the 2nd argument.+-- Allows for easy substring extraction, and retains the system of prefix/infix/suffix.+--+-- It is allowed to hand locations that only partially (or not at all) correspond to the @PIS@, but+-- then the resulting @PIS@ will be empty!++subPisLocation :: PIS i FwdPosition (BioSequence w) -> (FwdPosition, Int) -> PIS i FwdPosition (BioSequence w)+{-# Inline subPisLocation #-}+subPisLocation pis loc =+ let f z = subLocation z loc+ in over (pisPrefix._Just) f . over pisInfix f $ over (pisSuffix._Just) f pis++instance (Reversing (Location i FwdPosition (BioSequence w))) => Reversing (PIS i FwdPosition (BioSequence w)) where+ {-# Inline reversing #-}+ reversing pis+ = over (pisPrefix._Just) reversing . over pisInfix reversing . over (pisSuffix._Just) reversing+ . set pisPrefix (pis^.pisSuffix) . set pisSuffix (pis^.pisPrefix) $ pis++instance Complement (BioSequence w) => Complement (PIS i FwdPosition (BioSequence w)) where+ {-# Inline complement #-}+ complement =+ let f = over pisInfix (view complement) . over (pisPrefix._Just) (view complement) . over (pisSuffix._Just) (view complement)+ in iso f f++pisSequence :: Lens (PIS i p (BioSequence s)) (PIS i p (BioSequence t)) (BioSequence s) (BioSequence t)+{-# Inline pisSequence #-}+pisSequence = lens f t where+ v = view (locSequence.bioSequence)+ f (PIS p i s) = BioSequence $ maybe BS.empty v p `BS.append` v i `BS.append` maybe BS.empty v s+ t (PIS p i s) (BioSequence str) =+ let (pfx,ifxsfx) = over _1 BioSequence $ BS.splitAt (maybe 0 (BS.length . v) p) str+ (ifx,sfx ) = over both BioSequence $ BS.splitAt (BS.length $ v i) ifxsfx+ in PIS (set (_Just . locSequence) pfx p) (set locSequence ifx i) (set (_Just . locSequence) sfx s)++++-- | Given a @Location@ with a @BioSequence@, replace the sequence with its length.++locAsLength :: Location i FwdPosition (BioSequence w) -> Location i FwdPosition Int+{-# Inline locAsLength #-}+locAsLength = over locSequence (view (_BioSequence.to BS.length))++++-- | Provides a range in a notation as used by blast, for example. This+-- isomorphism can translate back as well. @FwdLocation - 8 4 ^. blastRange1 ==+-- 9 6 MinusStrand@, since these ranges are 1-based and start and end included.++blastRange1 :: (Location i FwdPosition Int) -> (Int,Int,Strand)+{-# Inline blastRange1 #-}+blastRange1 = f -- iso f t+ where+ f loc =+ let s = loc^.locPosition.fwdStart.to toInt1+ l = loc^.locSequence+ pm = loc^.locPosition.fwdStrand+ in case pm of PlusStrand -> (s,s+l,pm) ; MinusStrand -> (s+l,s,pm)+-- t (x,y,pm) =+-- let s = fromInt1 x+-- l = 1 + abs (x-y)+-- in Location (FwdPosition pm s) l++++-- | For each element, attach the prefix as well. The @Int@ indicates the maximal prefix length to+-- attach.+--+-- @1 2 3 4@ -> @01 12 23 34@+--+-- TODO are we sure this is correct for @MinusStrand@?++attachPrefixes+ :: ( Monad m, ModifyLocation p s )+ => Int+ -> SP.Stream (SP.Of (PIS i p s)) m r+ -> SP.Stream (SP.Of (PIS i p s)) m r+{-# Inlinable attachPrefixes #-}+attachPrefixes k = SP.map (\(Just w) -> w) . SP.drop 1 . SP.scan go Nothing id+ where+ go Nothing = Just+ go (Just p) = Just . set pisPrefix (Just . locTakeEnd k $ view pisInfix p)++++-- | For each element, attach the suffix as well.+--+-- @1 2 3 4@ -> @12 23 34 40@++attachSuffixes+ :: ( Monad m, ModifyLocation p s )+ => Int+ -> SP.Stream (SP.Of (PIS i p s)) m r+ -> SP.Stream (SP.Of (PIS i p s)) m r+{-# Inlinable attachSuffixes #-}+attachSuffixes k = loop Nothing+ where+ loop Nothing = \case+ SI.Return r -> SI.Return r+ SI.Effect m -> SI.Effect $ fmap (loop Nothing) m+ SI.Step (a SP.:> rest) -> loop (Just a) rest+ loop (Just p) = \case+ SI.Return r -> SI.Step (p SP.:> SI.Return r)+ SI.Effect m -> SI.Effect $ fmap (loop (Just p)) m+ SI.Step (a SP.:> rest) ->+ let p' = p & set pisSuffix (Just . locTake k $ view pisInfix a)+ in SI.Step (p' SP.:> loop (Just a) rest)+
− Biobase/Types/NumericalExtremes.hs
@@ -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 #-}-
− Biobase/Types/Odds.hs
@@ -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 #-}-
+ Biobase/Types/Position.hs view
@@ -0,0 +1,163 @@++-- | Annotate the genomic @position@ of features or elements. A @position@ has strand information,+-- and different ways to encode where a feature is located. The @position@ points to the first+-- element (e.g. nucleotide).+--+-- Together with the 'Biobase.Types.Location' module, it becomes possible to annotate substrings.++module Biobase.Types.Position where++import Control.DeepSeq+import Control.Lens hiding (Index, index)+import Data.Data+import GHC.Generics (Generic)+import GHC.TypeNats+import Prelude hiding (length)+import Text.Printf++import Biobase.Types.Index+import Biobase.Types.Strand+import Data.Info++{-++-- | Location information.++data Location = Location+ { _lStrand :: !Strand+ -- ^ On which strand are we+ , _lStart :: !(Index 0)+ -- ^ Start, 0-based+ , _lLength :: !Int+ -- ^ number of characters in this location+ , _lTotalLength :: !Int+ -- ^ the total length of the "contig" (or whatever) this location is positioned in.+ } deriving (Eq,Ord,Read,Show,Generic)+makeLenses ''Location+makePrisms ''Location++instance NFData Location++instance Semigroup Location where+ x <> y = let f z = z { _lLength = _lLength x + _lLength y }+ in case x^.lStrand of+ MinusStrand -> f y+ _otherStrand -> f x+ {-# Inline (<>) #-}++--instance Reversing Location where+-- {-# Inline reversing #-}+-- reversing = undefined+++-- | An isomorphism between locations, and triples of @Strand,Start,End@, where+-- end is inclusive. For @length==0@ locations, this will mean @start<end@ on+-- the plus strand.+--+-- This should hold for all @k@, in @Index k@.++startEndInclusive :: (KnownNat k) => Iso' Location (Strand, (Index k, Index k), Int)+{-# Inline startEndInclusive #-}+startEndInclusive = iso l2r r2l+ where l2r z = let s = z^.lStrand; f = z^.lStart; l = z^.lLength+ in (s, (reIndex f, reIndex $ f +. l -. 1), z^.lTotalLength)+ r2l (s,(f,t),ttl) = Location s (reIndex f) (delta f t + 1) ttl++-}++++-- | During streaming construction, it is possible that we know a feature is on the @-@ strand, but+-- the length of the contig is not known yet. In that case, 'FwdPosition' allows expressing the hit+-- in the coordinate system of the plus strand. Tools like blast do something similar, and express+-- locations on the minus as @y-x@ with @y>x@.+--+-- @+-- 0123456789+-- >-->+-- <--<+-- 9876543210+-- @+--+-- ++data FwdPosition+ -- | "Plus"-based location.+ = FwdPosition+ { _fwdStrand :: !Strand+ -- ^ Strand we are on+ , _fwdStart :: !(Index 0)+ -- ^ Start of the hit on the plus strand+ }+ deriving (Eq,Ord,Read,Show,Data,Typeable,Generic)+makeLenses ''FwdPosition+makePrisms ''FwdPosition++instance NFData FwdPosition++instance Info FwdPosition where+ info (FwdPosition s x) = printf "%s %d" (show s) (toInt0 x)++-- | Reversing a reversible location means moving the start to the end.++instance Reversing FwdPosition where+ {-# Inline reversing #-}+ reversing x = case x^.fwdStrand of+ PlusStrand -> set fwdStrand MinusStrand $ x+ MinusStrand -> set fwdStrand PlusStrand $ x+ UnknownStrand -> x++{-+++-- | Combining two FwdLocations yields the sum of their lengths. This assumes+-- that @x@ and @y@ are next to each other, or that it is ok if the @y@+-- @fwdStart@ information may be lost.+--+-- TODO provide associativity test in @properties@.++instance Semigroup FwdLocation where+ x <> y = over fwdLength (+ view fwdLength y) x+ {-# Inline (<>) #-}++instance ModifyLocation FwdLocation where+ locMoveLeftEnd k = over fwdStart (+. k) . over fwdLength (subtract k)+ locMoveRightEnd k = over fwdLength (+k)++-- | Given a location, take at most @k@ elements, and return a location after+-- this change.++fwdLocationTake :: Int -> FwdLocation -> FwdLocation+{-# Inline fwdLocationTake #-}+fwdLocationTake k' x =+ let l = x^.fwdLength+ k = max 0 $ min k' l -- deal with at most the length of the location+ in case x^.fwdStrand of+ MinusStrand -> set fwdLength k $ over fwdStart (+. (l-k)) x+ _otherStrand -> set fwdLength k $ x++-- | Given a location, drop at most @k@ elements, and return a location after+-- this change.+--+-- Note that @fwdLocationDrop 4 (FwdLocation PlusStrand 0 4) == FwdLocation 4 0@++fwdLocationDrop :: Int -> FwdLocation -> FwdLocation+{-# Inline fwdLocationDrop #-}+fwdLocationDrop k' x =+ let l = x^.fwdLength+ k = max 0 $ min k' l+ in case x^.fwdStrand of+ MinusStrand -> set fwdLength (l-k) $ x+ _otherStrand -> set fwdLength (l-k) $ over fwdStart (+. min k l) x++-- -- An isomorphism between a 'Location' and the pair @('FwdLocation',Int)@+-- -- exists.+-- +-- locationPartial :: Iso' Location (FwdLocation,Int)+-- {-# Inline locationPartial #-}+-- locationPartial = iso l2r r2l where+-- l2r l = undefined+-- r2l (p,z) = undefined++-}+
+ Biobase/Types/ReadingFrame.hs view
@@ -0,0 +1,54 @@++-- | Stranded reading frames.++module Biobase.Types.ReadingFrame where++import Control.Lens hiding (Index)+import GHC.Generics hiding (from)++import Biobase.Types.Index (Index, toInt0)+import Biobase.Types.Strand++++-- | 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.++newtype ReadingFrame = ReadingFrame { getReadingFrame ∷ Int }+ deriving (Eq,Ord,Generic,Show)+makeWrapped ''ReadingFrame++-- | Convert between @+1 ... +3@ and @ReadingFrame@.++rf ∷ Prism' Int ReadingFrame+{-# Inline rf #-}+rf = prism' getReadingFrame $ \k → let ak = abs k in+ if (ak <= 3 && ak >= 1) then Just (ReadingFrame k) else Nothing++-- | A lens for the strand++strandRF ∷ Lens' ReadingFrame Strand+{-# Inline strandRF #-}+strandRF = lens (\(ReadingFrame k) → if k < 0 then MinusStrand else PlusStrand)+ (\(ReadingFrame k) s → ReadingFrame $ if s == PlusStrand then abs k else (negate $ abs k))++-- |+--+-- @pred@ and @succ@ are correct, if the input is a legal 'ReadingFrame'.++instance Enum ReadingFrame where+ {-# Inline toEnum #-}+ toEnum k = case k^?rf of Just rf → rf ; Nothing → error $ show k ++ " is not a legal reading frame"+ {-# Inline fromEnum #-}+ fromEnum = getReadingFrame++-- |+--+-- 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+
− Biobase/Types/Sequence.hs
@@ -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.--- <http://hyperphysics.phy-astr.gsu.edu/hbase/Organic/transcription.html>------ @@ 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 #-}-
+ Biobase/Types/Shape.hs view
@@ -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+
Biobase/Types/Strand.hs view
@@ -1,63 +1,93 @@ -- | Strand information. A newtyped version, complete with serialization,--- pattern synonyms, being a @PrimitiveArray@ index type, etc.+-- pattern synonyms, being a @PrimitiveArray@ index type, etc. The strand+-- information includes @+@, @-@, as well as the (GFF3) @.@ not stranded, and+-- @?@ for unknown strand information. ----- TODO will be expanded to encode biological sense information more--- clearly: <http://en.wikipedia.org/wiki/Sense_%28molecular_biology%29>.+-- TODO will be expanded to encode biological sense information more clearly:+-- <http://en.wikipedia.org/wiki/Sense_%28molecular_biology%29>. module Biobase.Types.Strand where import Control.DeepSeq+import Control.Lens hiding (Index) import Control.Monad (guard) import Data.Aeson import Data.Binary+import Data.Data (Data) import Data.Hashable (Hashable) import Data.Serialize (Serialize)-import Data.Vector.Fusion.Stream.Monadic (Step(..))+import Data.Typeable (Typeable)+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)+ deriving (Eq,Ord,Generic,Data,Typeable) instance Show Strand where- show P = "+"- show M = "-"+ show PlusStrand = "PlusStrand"+ show MinusStrand = "MinusStrand"+ show NotStranded = "NotStranded"+ show UnknownStrand = "UnknownStrand" 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)+ "NotStranded" -> return (NotStranded, s)+ "UnknownStrand" -> return (UnknownStrand, s)+ [x] | x `elem` ("+Pp" :: String) -> return (PlusStrand,s)+ | x `elem` ("-Mm" :: String) -> return (MinusStrand,s)+ | x `elem` ("." :: String) -> return (NotStranded,s)+ | x `elem` ("?" :: String) -> return (UnknownStrand,s)+ _ -> [] instance Bounded Strand where- minBound = P- maxBound = M+ minBound = PlusStrand+ maxBound = UnknownStrand instance Enum Strand where- succ P = M- succ M = error "succ M"- pred M = P- pred P = error "pred P"- toEnum i | i>=0 && i<=1 = Strand i+ succ (Strand k)+ | k < 0 = error "succ undefined strand"+ | k == 3 = error "succ UnknownStrand"+ | k > 3 = error "succ undefined strand"+ | otherwise = Strand (k+1)+ pred (Strand k)+ | k < 0 = error "pred undefined strand"+ | k == 0 = error "pred UnknownStrand"+ | k > 3 = error "pred undefined strand"+ | otherwise = Strand (k-1)+ toEnum i | i>=0 && i<=3 = Strand i toEnum i = error $ "toEnum (Strand)" ++ show i fromEnum = getStrand -pattern P = Strand 0-pattern M = Strand 1+instance Reversing Strand where+ reversing PlusStrand = MinusStrand+ reversing MinusStrand = PlusStrand+ reversing x = x -pattern Sense = P-pattern AntiSense = M+pattern PlusStrand = Strand 0+pattern MinusStrand = Strand 1+pattern NotStranded = Strand 2+pattern UnknownStrand = Strand 3 +-- TODO Sense and Antisense are somewhat different++--pattern Sense = P+--pattern AntiSense = M+ instance Binary Strand instance Serialize Strand instance ToJSON Strand@@ -69,19 +99,26 @@ [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 #-}+ fromLinearIndex _ = Strand+ {-# Inline [0] fromLinearIndex #-}+ showBound (LtStrand k) = ["LtStrand " ++ show k]+ showIndex (Strand k) = ["Strand " ++ show k] 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 +126,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@@ -98,11 +135,11 @@ {-# Inline [0] step #-} {-# Inline streamDown #-} -instance IndexStream Strand+-- instance IndexStream Strand instance Arbitrary Strand where arbitrary = do- b <- choose (0,1)+ b <- choose (0,3) return $ Strand b shrink (Strand j) | 0<j = [Strand $ j-1]
Biobase/Types/Structure.hs view
@@ -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
Biobase/Types/Taxonomy.hs view
@@ -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
BiobaseTypes.cabal view
@@ -1,17 +1,17 @@+cabal-version: 2.2 name: BiobaseTypes-version: 0.1.3.0-author: Christian Hoener zu Siederdissen, 2015 - 2017-copyright: Christian Hoener zu Siederdissen, 2015 - 2017+version: 0.2.1.0+author: Christian Hoener zu Siederdissen, 2015 - 2021+copyright: Christian Hoener zu Siederdissen, 2015 - 2021 homepage: https://github.com/choener/BiobaseTypes bug-reports: https://github.com/choener/BiobaseTypes/issues maintainer: choener@bioinf.uni-leipzig.de category: Data Structures, Bioinformatics-license: BSD3+license: BSD-3-Clause license-file: LICENSE build-type: Simple stability: experimental-cabal-version: >= 1.10.0-tested-with: GHC == 8.0.2, GHC == 8.2.1+tested-with: GHC == 8.8, GHC == 8.10, GHC == 9.0 synopsis: Collection of types for bioinformatics description: Types used in a number of bioinformatics libraries.@@ -20,8 +20,6 @@ . * energies .- * numerics- . * biostring wrappers @@ -32,9 +30,10 @@ -library+common deps build-depends: base >= 4.7 && < 5.0 , aeson >= 0.8+ , attoparsec >= 0.13 , binary >= 0.7 , bytestring , cereal >= 0.4@@ -49,6 +48,7 @@ , mtl , primitive >= 0.5 , QuickCheck >= 2.7+ , streaming >= 0.1 , string-conversions >= 0.4 , text >= 1.0 , text-binary >= 0.2@@ -58,28 +58,19 @@ , vector-th-unbox >= 0.2 -- , bimaps == 0.1.0.*- , PrimitiveArray == 0.8.0.*- exposed-modules:- Biobase.Types.Accession- Biobase.Types.Bitscore- Biobase.Types.Energy- Biobase.Types.Evalue- Biobase.Types.Index- Biobase.Types.Index.Type- Biobase.Types.Names- Biobase.Types.Names.Internal- Biobase.Types.NumericalExtremes- Biobase.Types.Odds- Biobase.Types.Sequence- Biobase.Types.Strand- Biobase.Types.Structure- Biobase.Types.Taxonomy+ , DPutils == 0.1.1.*+ , ForestStructures == 0.0.1.*+ , PrimitiveArray >= 0.10.1.1 && < 0.10.2+ , SciBaseTypes == 0.1.1.* default-language: Haskell2010 default-extensions: BangPatterns , DataKinds , DeriveDataTypeable+ , DeriveFoldable , DeriveGeneric+ , DeriveTraversable+ , DerivingStrategies , FlexibleContexts , FlexibleInstances , GeneralizedNewtypeDeriving@@ -88,33 +79,64 @@ , MultiParamTypeClasses , OverloadedStrings , PatternSynonyms+ , PolyKinds+ , RankNTypes+ , RecordWildCards , ScopedTypeVariables , StandaloneDeriving , TemplateHaskell+ , TypeApplications , TypeFamilies , TypeOperators+ , TupleSections+ , UndecidableInstances , UnicodeSyntax+ , ViewPatterns ghc-options: -O2 -funbox-strict-fields +library+ import:+ deps+ exposed-modules:+ Biobase.Types.Accession+ Biobase.Types.BioSequence+ Biobase.Types.Bitscore+ Biobase.Types.Codon+ Biobase.Types.Energy+ Biobase.Types.Evalue+ Biobase.Types.Index+ Biobase.Types.Index.Type+ Biobase.Types.Location+ Biobase.Types.Names+ Biobase.Types.Names.Internal+ Biobase.Types.Position+ Biobase.Types.ReadingFrame+ Biobase.Types.Shape+ Biobase.Types.Strand+ Biobase.Types.Structure+ Biobase.Types.Taxonomy+ DP.Backtraced.BioSequence+ DP.Backtraced.Codon+++ test-suite properties+ import:+ deps type: exitcode-stdio-1.0 main-is: properties.hs- ghc-options:- -threaded -rtsopts -with-rtsopts=-N+-- ghc-options:+-- -threaded -rtsopts -with-rtsopts=-N hs-source-dirs: tests- default-language:- Haskell2010- default-extensions: ScopedTypeVariables- , TemplateHaskell build-depends: base- , QuickCheck , tasty >= 0.11+ , tasty-hunit >= 0.10 , tasty-quickcheck >= 0.8 , tasty-th >= 0.1 --
+ DP/Backtraced/BioSequence.hs view
@@ -0,0 +1,3 @@++module DP.Backtraced.BioSequence where+
+ DP/Backtraced/Codon.hs view
@@ -0,0 +1,46 @@++-- | The Backtraced column structure is for codon-based alignments, including+-- special cases.++module DP.Backtraced.Codon where++import Data.ByteString (ByteString)+import Data.Vector (Vector)+import GHC.Generics (Generic)++import Biobase.Types.Codon++++-- | A single 'Backtraced' column. Since such a column will be part of a+-- @Backtraced (Z:.BtCodon c aa:. ...)@ structure, it is always possible to+-- extend even further, by having more entries.++data BtCodon c aa+ -- | A canonical match. A codon and the translated amino acid need to be set.+ = Match+ { _codon ∷ !(Codon c)+ , _aa ∷ !aa+ }+ -- | A frameshifting match. The vector of frameshifted nucleotides will have+ -- a number of characters @c@, that encode for a single amino acid.+ | Frameshift+ { _frameshift ∷ !(Vector c)+ , _aa ∷ !aa+ }+ | Insert+ { _codon ∷ !(Codon c)+ , _aa ∷ !aa+ }+ | Shifted+ { _frameshift ∷ !(Vector c)+ , _aa ∷ !aa+ }+ | Region+ { _region ∷ !(Vector c)+ , _annotation ∷ !ByteString+ }+ | Delete+ {+ }+
LICENSE view
@@ -1,4 +1,4 @@-Copyright Christian Hoener zu Siederdissen 2015+Copyright Christian Hoener zu Siederdissen 2015-2019 All rights reserved.
README.md view
@@ -1,4 +1,5 @@-[](https://travis-ci.org/choener/BiobaseTypes)++ # BiobaseTypes
changelog.md view
@@ -1,3 +1,25 @@+0.2.1.0+-------++- CI/hackage github actions+- dependency updates++0.2.0.1+-------++- minor version bumped due to OrderedBits++0.2.0.0+-------++- unified treatment of bio sequences with one phantom-typed newtype.++0.1.4.0+-------++- changes to indexing and others+- some changes that probably require a bump+ 0.1.3.0 -------
tests/properties.hs view
@@ -1,16 +1,25 @@ 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.QuickCheck.Property () import Test.Tasty+import Test.Tasty.HUnit import Test.Tasty.QuickCheck (testProperty) import Test.Tasty.TH +import Biobase.Types.BioSequence import Biobase.Types.Bitscore-import Biobase.Types.NumericalExtremes-+import Biobase.Types.Location+import Biobase.Types.Shape+import Biobase.Types.Strand+import Biobase.Types.Structure+import Biobase.Types.Index as I +{- -- * Bitscore conversions @@ -19,8 +28,107 @@ --prop_ScoreProb (Positive null) x = Bitscore x ~= prob2Score null (score2Prob null $ Bitscore x) ++-- * sequence properties++-- complement twice++prop_complement_twice_DNA (dna ∷ BioSequence DNA) = dna == dna^.complement.complement++prop_complement_twice_RNA (rna ∷ BioSequence RNA) = rna == rna^.complement.complement++prop_transcribe_twice_DNA (dna ∷ BioSequence DNA) = 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++prop_FwdLocationPlusTake (NonNegative (p ∷ Int), NonNegative (l ∷ Int), NonNegative (k ∷ Int))+ | check = True+ | otherwise = traceShow (p,l,k,fwdloc,taken,manual) check+ where fwdloc = FwdLocation PlusStrand (I.index p) l+ check = taken == manual+ taken = fwdLocationTake k fwdloc+ manual = FwdLocation PlusStrand (I.index p) (max 0 $ min l k)++prop_FwdLocationPlusDrop (NonNegative (p ∷ Int), NonNegative (l ∷ Int), NonNegative (k ∷ Int))+ | check = True+ | otherwise = traceShow (p,l,k,fwdloc,dropped,manual) check+ where fwdloc = FwdLocation PlusStrand (I.index p) l+ check = dropped == manual+ dropped = fwdLocationDrop k fwdloc+ manual = FwdLocation PlusStrand (I.index $ p + min l k) (max 0 $ l-k)++-- | Given a BioSequenceWindow, and different takes and drops, check wether what we have corresponds to what we want++case_bswTakeDrop ∷ Assertion+case_bswTakeDrop = do+ let wp = BioSequenceWindow @"DNA" @DNA "test" 1 "ACGTAC" 3 (FwdLocation PlusStrand 0 6)+ wm = BioSequenceWindow @"DNA" @DNA "test" 3 "CATGCA" 1 (FwdLocation MinusStrand 0 6)+ --+ bswTake 0 wp @?= BioSequenceWindow "test" 0 "" 0 (FwdLocation PlusStrand 0 0)+ bswTake 1 wp @?= BioSequenceWindow "test" 1 "A" 0 (FwdLocation PlusStrand 0 1)+ bswTake 2 wp @?= BioSequenceWindow "test" 1 "AC" 0 (FwdLocation PlusStrand 0 2)+ bswTake 6 wp @?= BioSequenceWindow "test" 1 "ACGTAC" 3 (FwdLocation PlusStrand 0 6)+ --+ bswDrop 0 wp @?= BioSequenceWindow "test" 1 "ACGTAC" 3 (FwdLocation PlusStrand 0 6)+ bswDrop 1 wp @?= BioSequenceWindow "test" 0 "CGTAC" 3 (FwdLocation PlusStrand 1 5)+ bswDrop 6 wp @?= BioSequenceWindow "test" 0 "" 0 (FwdLocation PlusStrand 6 0)+ --+ bswTake 0 wm @?= BioSequenceWindow "test" 0 "" 0 (FwdLocation MinusStrand 6 0)+ bswTake 1 wm @?= BioSequenceWindow "test" 1 "C" 0 (FwdLocation MinusStrand 5 1)+ bswTake 2 wm @?= BioSequenceWindow "test" 2 "CA" 0 (FwdLocation MinusStrand 4 2)+ bswTake 3 wm @?= BioSequenceWindow "test" 3 "CAT" 0 (FwdLocation MinusStrand 3 3)+ bswTake 4 wm @?= BioSequenceWindow "test" 3 "CATG" 0 (FwdLocation MinusStrand 2 4)+ bswTake 5 wm @?= BioSequenceWindow "test" 3 "CATGC" 0 (FwdLocation MinusStrand 1 5)+ bswTake 6 wm @?= BioSequenceWindow "test" 3 "CATGCA" 1 (FwdLocation MinusStrand 0 6)+ --+ bswDrop 0 wm @?= BioSequenceWindow "test" 3 "CATGCA" 1 (FwdLocation MinusStrand 0 6)+ bswDrop 1 wm @?= BioSequenceWindow "test" 2 "ATGCA" 1 (FwdLocation MinusStrand 0 5)+ bswDrop 2 wm @?= BioSequenceWindow "test" 1 "TGCA" 1 (FwdLocation MinusStrand 0 4)+ bswDrop 5 wm @?= BioSequenceWindow "test" 0 "A" 1 (FwdLocation MinusStrand 0 1)+ bswDrop 6 wm @?= BioSequenceWindow "test" 0 "" 0 (FwdLocation MinusStrand 0 0)+ --+ -- TODO consider having [take,take,drop,drop], generate all permutations;+ -- they should all yield the same result.+ --+++++-- * generic stuff+ a ~= b = abs (b-a) <= 10e-6 main :: IO () main = $(defaultMainGenerator)++-}++main :: IO ()+main = return ()