diff --git a/Biobase/Primary.hs b/Biobase/Primary.hs --- a/Biobase/Primary.hs +++ b/Biobase/Primary.hs @@ -1,219 +1,23 @@ -{-# LANGUAGE PackageImports #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE OverlappingInstances #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE StandaloneDeriving #-} --- | The primary structure: interface to efficient encoding of RNA and DNA --- sequences. The design aims toward the 'vector' library and repa. In --- particular, everything is strict; if you want to stream full genomes, use --- 'text' or lazy 'bytestring's instead and cast to Biobase.Primary definitions --- only at the last moment. +-- | -- --- NOTE individual nucleotides are encoded is 'Int's internally without any --- tagging. This means that we have no way of deciding if we are dealing with --- RNA or DNA on this level. - -module Biobase.Primary where - -import "PrimitiveArray" Data.Array.Repa.Index -import "PrimitiveArray" Data.Array.Repa.Shape -import Data.Char (toUpper) -import Data.ExtShape -import Data.Ix (Ix(..)) -import Data.Primitive.Types -import Data.Tuple (swap) -import GHC.Base (remInt,quotInt) -import qualified Data.ByteString.Char8 as BS -import qualified Data.ByteString.Lazy.Char8 as BSL -import qualified Data.Text as T -import qualified Data.Vector.Generic as VG -import qualified Data.Vector.Generic.Mutable as VGM -import qualified Data.Vector.Unboxed as VU - -import Data.PrimitiveArray -import Data.PrimitiveArray.Unboxed.Zero - -import Biobase.Primary.Bounds - - - --- * Convert different types of sequence representations to the internal --- "Primary Structure" representation - --- | Given a sequence of nucleotides encoded in some "text-form", create a --- 'Nuc'-based unboxed vector. - -class MkPrimary a where - mkPrimary :: a -> Primary - -type Primary = Arr0 DIM1 Nuc - -instance Eq Primary where - xs == ys - | bx==by = sliceEq xs zeroDim ys zeroDim bx - | otherwise = False - where (_,bx) = bounds xs - (_,by) = bounds ys - -instance Ord Primary where - xs <= ys - | bx==by = toList xs <= toList ys - | otherwise = bx<=by - where (_,Z:.bx) = bounds xs - (_,Z:.by) = bounds ys - - - --- * Efficient nucleotide encoding - --- A 'Nuc'leotide is simply an Int wrapped up. 'nIMI' provides for --- intermolecular initialization, 'nN' stands for "any" nucleotides, 'nA', --- 'nC', 'nG', 'nT' / 'nU' are normal nucleotides. - -newtype Nuc = Nuc {unNuc :: Int} - deriving (Eq,Ord,Ix) - -(nN : nA : nC : nG : nT : nIMI : nUndefined : _) = map Nuc [0 .. ] -nU = nT - -acgt = [nA..nT] -acgu = acgt -nacgt = [nN..nT] -nacgu = nacgt - --- | Translate between 'Char's and 'Nuc's. - -mkNuc :: Char -> Nuc -mkNuc = f . toUpper where - f k - | Just v <- k `lookup` charNucList = v - | otherwise = nN - -fromNuc :: Nuc -> Char -fromNuc = f where - f k - | Just v <- k `lookup` nucCharList = v - | otherwise = 'N' - -charNucList = - [ ('N',nN) - , ('A',nA) - , ('C',nC) - , ('G',nG) - , ('T',nT) - , ('U',nU) - ] - -nucCharList = map swap charNucList - - - --- * Instances of different type classes - --- ** instances for 'Nuc' - --- | Human-readable Show instance. - -instance Show Nuc where - show n = [fromNuc n] - --- | Human-readable Read instance. - -instance Read Nuc where - readsPrec p [] = [] - readsPrec p (x:xs) - | x ==' ' = readsPrec p xs - | Just n <- x `lookup` charNucList = [(n,xs)] - | otherwise = [] - --- for vectors - -deriving instance Prim Nuc -deriving instance VGM.MVector VU.MVector Nuc -deriving instance VG.Vector VU.Vector Nuc -deriving instance VU.Unbox Nuc - --- Shape-based indexing. Nucleotide representations go from nN (0) to nU (4), --- with additional symbols being available for specialized problems. This is a --- bit of a problem for shape-based indexing. In particular, we need to be --- careful with size operations. To include, say, all of nN to nU one needs a --- size of (z:.nIMI), as nIMI is the first element not in the size anymore. - -instance (Shape sh,Show sh) => Shape (sh :. Nuc) where - rank (sh:._) = rank sh + 1 - zeroDim = zeroDim:.Nuc 0 - unitDim = unitDim:.Nuc 1 -- TODO does this one make sense? - intersectDim (sh1:.n1) (sh2:.n2) = intersectDim sh1 sh2 :. min n1 n2 - addDim (sh1:.Nuc n1) (sh2:.Nuc n2) = addDim sh1 sh2 :. Nuc (n1+n2) -- TODO will not necessarily yield a valid Nuc - size (sh1:.Nuc n) = size sh1 * n - sizeIsValid (sh1:.Nuc n) = sizeIsValid (sh1:.n) - toIndex (sh1:.Nuc sh2) (sh1':.Nuc sh2') = toIndex (sh1:.sh2) (sh1':.sh2') - fromIndex (ds:.Nuc d) n = fromIndex ds (n `quotInt` d) :. Nuc r where - r | rank ds == 0 = n - | otherwise = n `remInt` d - inShapeRange (sh1:.n1) (sh2:.n2) (idx:.i) = i>=n1 && i error "empty list in shapeOfList/Primary" - x:xs -> shapeOfList xs :. Nuc x - deepSeq (sh:.n) x = deepSeq sh (n `seq` x) - {-# INLINE rank #-} - {-# INLINE zeroDim #-} - {-# INLINE unitDim #-} - {-# INLINE intersectDim #-} - {-# INLINE addDim #-} - {-# INLINE size #-} - {-# INLINE sizeIsValid #-} - {-# INLINE toIndex #-} - {-# INLINE fromIndex #-} - {-# INLINE inShapeRange #-} - {-# INLINE listOfShape #-} - {-# INLINE shapeOfList #-} - {-# INLINE deepSeq #-} - -instance (Shape sh, Show sh, ExtShape sh) => ExtShape (sh :. Nuc) where - subDim (sh1:.Nuc n1) (sh2:.Nuc n2) = subDim sh1 sh2 :. Nuc (n1-n2) - rangeList (sh1:.Nuc n1) (sh2:.Nuc n2) = [ sh:.Nuc n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2)]] - --- | The bounded instance from GHC proper. Captures all defined symbols. - -instance Bounded Nuc where - minBound = nN - maxBound = nT - --- | Special bounds for energy / score arrays - -instance Bounds Nuc where - minNormal = nA - maxNormal = nT - minExtended = nN - maxExtended = nT - --- | Enum - -instance Enum Nuc where - toEnum = Nuc - fromEnum = unNuc - --- ** Instances for 'MkPrimary' - -instance MkPrimary String where - mkPrimary xs = fromList (Z:.0) (Z:.length xs -1) $ map mkNuc xs - -instance MkPrimary BS.ByteString where - mkPrimary = mkPrimary . BS.unpack - -instance MkPrimary BSL.ByteString where - mkPrimary = mkPrimary . BSL.unpack +-- TODO make sequence types 'stringable'? -instance MkPrimary T.Text where - mkPrimary = mkPrimary . T.unpack +module Biobase.Primary + ( module Biobase.Primary.AA + , module Biobase.Primary.Hashed + , module Biobase.Primary.IUPAC + , module Biobase.Primary.Letter + , module Biobase.Primary.Nuc + , module Biobase.Primary.Trans + , module Biobase.Primary.Unknown + ) where -instance MkPrimary [Nuc] where - mkPrimary xs = fromList (Z:.0) (Z:.length xs -1) xs +import Biobase.Primary.AA hiding (Stop,A,B,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,X,Y,Z,Undef) +import Biobase.Primary.Hashed +import Biobase.Primary.IUPAC hiding (A,C,G,T,U,W,S,M,K,R,Y,B,D,H,V,N) +import Biobase.Primary.Letter +import Biobase.Primary.Nuc +import Biobase.Primary.Trans +import Biobase.Primary.Unknown diff --git a/Biobase/Primary/AA.hs b/Biobase/Primary/AA.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/AA.hs @@ -0,0 +1,171 @@ + +-- | This module has the translation tables for the genetic code. +-- +-- In addition, @Any@ is included to denote that any amino acid is ok, and +-- @Unknown@ to denote unknown data. We do have a symbol 'Undef' for undefined +-- amino acids, which denotes error condition. +-- +-- TODO this nomenclature might change! + +module Biobase.Primary.AA where + +import Control.Arrow ((***),first) +import Data.Aeson +import Data.Hashable +import Data.Ix (Ix(..)) +import Data.Map.Strict (Map) +import Data.Primitive.Types +import Data.Tuple (swap) +import Data.Vector.Unboxed.Deriving +import GHC.Base (remInt,quotInt) +import GHC.Generics (Generic) +import qualified Data.Bijection.HashMap as B +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as BSL +import qualified Data.Foldable as F +import qualified Data.Map.Strict as M +import qualified Data.Text as T +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Unboxed as VU +import qualified GHC.Exts as GHC + +import Biobase.Types.BioSequence +import Data.Info + +import Biobase.Primary.Letter + + + +pattern Stop = Letter 0 :: Letter AA n +pattern A = Letter 1 :: Letter AA n +pattern B = Letter 2 :: Letter AA n +pattern C = Letter 3 :: Letter AA n +pattern D = Letter 4 :: Letter AA n +pattern E = Letter 5 :: Letter AA n +pattern F = Letter 6 :: Letter AA n +pattern G = Letter 7 :: Letter AA n +pattern H = Letter 8 :: Letter AA n +pattern I = Letter 9 :: Letter AA n +pattern K = Letter 10 :: Letter AA n +pattern L = Letter 11 :: Letter AA n +pattern M = Letter 12 :: Letter AA n +pattern N = Letter 13 :: Letter AA n +pattern P = Letter 14 :: Letter AA n +pattern Q = Letter 15 :: Letter AA n +pattern R = Letter 16 :: Letter AA n +pattern S = Letter 17 :: Letter AA n +pattern T = Letter 18 :: Letter AA n +pattern V = Letter 19 :: Letter AA n +pattern W = Letter 20 :: Letter AA n +pattern X = Letter 21 :: Letter AA n +pattern Y = Letter 22 :: Letter AA n +pattern Z = Letter 23 :: Letter AA n +pattern Any = Letter 24 :: Letter AA n -- TODO @Any == X@ supposedly! +pattern Unknown = Letter 25 :: Letter AA n +pattern Undef = Letter 26 :: Letter AA n + +-- * Creating functions and aa data. + +aa :: Int -> Letter AA n +aa = Letter +{-# Inline aa #-} + +aaRange = VU.fromList [Stop .. pred Undef] +{-# NoInline aaRange #-} + +instance Bounded (Letter AA n) where + minBound = Stop + maxBound = Undef + +instance LetterChar AA n where + letterChar = aaChar + charLetter = charAA + +instance ToJSON (Letter AA n) where + toJSON = toJSON . letterChar + +instance FromJSON (Letter AA n) where + parseJSON = fmap charLetter . parseJSON + +instance Info (Letter AA n) where + info = (:[]) . aaChar + +--instance (GHC.IsString f) => ToJSON (Pretty f (Letter AA)) where +-- toJSON = toJSON . T.pack . map letterChar . GHC.toList . getPretty + +-- | Translate 'Char' amino acid representation into efficient 'AA' newtype. + +charAA :: Char -> Letter AA n +charAA = B.findWithDefaultL Undef charBaa +{-# INLINE charAA #-} + +-- | 'Char' representation of an 'AA'. + +aaChar :: Letter AA n -> Char +aaChar = B.findWithDefaultR '?' charBaa +{-# INLINE aaChar #-} + +-- * lookup tables + +charBaa :: B.Bimap (B.HashMap Char (Letter AA n)) (B.HashMap (Letter AA n) Char) +charBaa = B.fromList + [ ('*',Stop) + , ('A',A) + , ('B',B) + , ('C',C) + , ('D',D) + , ('E',E) + , ('F',F) + , ('G',G) + , ('H',H) + , ('I',I) + , ('K',K) + , ('L',L) + , ('M',M) + , ('N',N) + , ('P',P) + , ('Q',Q) + , ('R',R) + , ('S',S) + , ('T',T) + , ('V',V) + , ('W',W) + , ('X',X) + , ('Y',Y) + , ('Z',Z) + , ('?',Unknown) + ] +{-# NOINLINE charBaa #-} + +-- | List of the twenty "default" amino acids. Used, for example, by HMMer. + +twentyAA :: VU.Vector (Letter AA n) +twentyAA = VU.fromList [ A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y ] +{-# NoInline twentyAA #-} + + +-- * instances + +instance Show (Letter AA n) where + show n = [aaChar n] + +instance Read (Letter AA n) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | aa <- charAA x = [(aa,xs)] + | otherwise = [] + +instance Enum (Letter AA n) where + succ Undef = error "succ/Undef:AA" + succ (Letter x) = Letter $ x+1 + pred Stop = error "pred/Stop:AA" + pred (Letter x) = Letter $ x-1 + toEnum k | k>=0 && k<=(getLetter Undef) = Letter k + toEnum k = error $ "toEnum/Letter RNA " ++ show k + fromEnum (Letter k) = k + +instance MkPrimary (VU.Vector Char) AA n where + primary = VU.map charAA + diff --git a/Biobase/Primary/Bounds.hs b/Biobase/Primary/Bounds.hs --- a/Biobase/Primary/Bounds.hs +++ b/Biobase/Primary/Bounds.hs @@ -5,7 +5,7 @@ module Biobase.Primary.Bounds where - +{- -- | 'minNormal' and 'maxNormal' encode for, say, ACGU; while 'minExtended' and -- 'maxExtended' would allow 'N' as well. See Biobase.RNA and @@ -50,3 +50,6 @@ maxNormal = (maxNormal, maxNormal, maxNormal, maxNormal, maxNormal, maxNormal) minExtended = (minExtended, minExtended, minExtended, minExtended, minExtended, minExtended) maxExtended = (maxExtended, maxExtended, maxExtended, maxExtended, maxExtended, maxExtended) + +-} + diff --git a/Biobase/Primary/Hashed.hs b/Biobase/Primary/Hashed.hs --- a/Biobase/Primary/Hashed.hs +++ b/Biobase/Primary/Hashed.hs @@ -1,50 +1,46 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} --- | Fast hash functions for 'Primary' sequences. A hash is just an 'Int', so --- use these only for short sequences. +-- | Fast hash functions for 'Primary' sequences. This function maps +-- primary sequences to a continuous set of Ints @[0 ..]@ where the maximum +-- is dependent on the input length. This allows us to map short sequences +-- into contiguous memory locations. Useful for, say, energy lookup tables. module Biobase.Primary.Hashed where -import Control.Exception.Base (assert) -import Data.Ix -import Data.Primitive.Types +import Data.Ix +import Data.Primitive.Types +import Data.Vector.Unboxed.Deriving import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU -import Data.PrimitiveArray +import Biobase.Primary.Letter -import Biobase.Primary +-- | The hash of a primary sequence. -newtype HashedPrimary = HashedPrimary Int +newtype HashedPrimary t n = HashedPrimary { unHashedPrimary :: Int } deriving (Eq,Ord,Ix,Read,Show,Enum,Bounded) -deriving instance Prim HashedPrimary -deriving instance VGM.MVector VU.MVector HashedPrimary -deriving instance VG.Vector VU.Vector HashedPrimary -deriving instance VU.Unbox HashedPrimary +derivingUnbox "HashedPrimary" + [t| forall t n . HashedPrimary t n -> Int |] [| unHashedPrimary |] [| HashedPrimary |] -- | Given a piece of primary sequence information, reduce it to an index. --- --- Will throw an assertion in debug code if 'ps' are not within bounds. Note --- that "mkPrimary [minBound]" and "mkPrimary [minBound,minBound]" map to the --- same index. Meaning that indices are only unique within the same length --- group. Furthermore, indices with different (l,u)-bounds are not compatible --- with each other. All indices start at 0. --- -- The empty input produces an index of 0. --- --- TODO currently goes the very inefficient way of creating a temporary vector --- for 'ps'. We could in O(1) create a vector from a Primary ... -mkHashedPrimary :: (Nuc,Nuc) -> Primary -> HashedPrimary -mkHashedPrimary (l,u) ps' = assert (VU.all (\p -> l<=p && p<=u) ps) $ HashedPrimary idx where - idx = VU.sum $ VU.zipWith f ps (VU.enumFromStepN (VU.length ps -1) (-1) (VU.length ps)) - f p c = (unNuc p - unNuc l) * (cnst^c) - cnst = unNuc u - unNuc l + 1 - ps = VU.fromList $ toList ps' +mkHashedPrimary :: forall t n . (VU.Unbox (Letter t n), Bounded (Letter t n), Enum (Letter t n)) => Primary t n -> HashedPrimary t n +mkHashedPrimary = HashedPrimary . fst . VU.foldl' f (0, 1) where + f (z, c) n = (z + c * (fromEnum n +1), c * (fromEnum (maxBound :: Letter t n) + 1)) {-# INLINE mkHashedPrimary #-} + +-- | Turn a hash back into a sequence. Will fail if the resulting sequence +-- has more than 100 elements. + +hash2primary :: forall t n . (VU.Unbox (Letter t n), Bounded (Letter t n), Enum (Letter t n)) => HashedPrimary t n -> Primary t n +hash2primary (HashedPrimary h) = VU.unfoldrN l f h where + m = fromEnum (maxBound :: Letter t n) +1 + l = VU.length . VU.takeWhile (>0) . VU.iterateN 100 (`div` m) $ h + f k = if k>0 then Just (toEnum $ ((k-1) `mod` m) , (k-1) `div` m) + else Nothing +{-# INLINE hash2primary #-} + diff --git a/Biobase/Primary/IUPAC.hs b/Biobase/Primary/IUPAC.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/IUPAC.hs @@ -0,0 +1,155 @@ + +-- | Degenerate base symbol representation. We use the same conventions as in +-- <> which ignores +-- @U@racil, except if it stands alone for @Char@ and @XNA@ targets. If the +-- 'Degenerate' target is @RNA@, then we create @U@s instead of @T@s. +-- +-- TODO Shall we handle 'Complement' for degenerates? + +module Biobase.Primary.IUPAC where + +import Control.Arrow ((***)) +import Data.ByteString.Char8 (ByteString,unpack) +import Data.Char (toUpper) +import Data.FileEmbed (makeRelativeToProject, embedFile) +import Data.List (nub,sort) +import Data.String +import Data.Tuple (swap) +import qualified Data.Vector.Unboxed as VU +import Control.Category ((>>>)) + +import Biobase.Types.BioSequence + +import Biobase.Primary.Letter +import Biobase.Primary.Nuc +import qualified Biobase.Primary.Nuc.RNA as R + + +-- | Allow the full, including degenerates, alphabet. + +data DEG + +pattern A = Letter 0 :: Letter DEG n +pattern C = Letter 1 :: Letter DEG n +pattern G = Letter 2 :: Letter DEG n +pattern T = Letter 3 :: Letter DEG n +pattern U = Letter 4 :: Letter DEG n +pattern W = Letter 5 :: Letter DEG n +pattern S = Letter 6 :: Letter DEG n +pattern M = Letter 7 :: Letter DEG n +pattern K = Letter 8 :: Letter DEG n +pattern R = Letter 9 :: Letter DEG n +pattern Y = Letter 10 :: Letter DEG n +pattern B = Letter 11 :: Letter DEG n +pattern D = Letter 12 :: Letter DEG n +pattern H = Letter 13 :: Letter DEG n +pattern V = Letter 14 :: Letter DEG n +pattern N = Letter 15 :: Letter DEG n + +instance Bounded (Letter DEG n) where + minBound = A + maxBound = N + +instance Enum (Letter DEG n) where + succ N = error "succ/N:DEG" + succ (Letter x) = Letter $ x+1 + pred A = error "pred/A:DEG" + pred (Letter x) = Letter $ x-1 + toEnum k | k>=0 && k<=15 = Letter k + toEnum k = error $ "toEnum/Letter DEG " ++ show k + fromEnum (Letter k) = k + +charDEG = toUpper >>> \case + 'A' -> A + 'C' -> C + 'G' -> G + 'T' -> T + 'U' -> U + 'W' -> W + 'S' -> S + 'M' -> M + 'K' -> K + 'R' -> R + 'Y' -> Y + 'B' -> B + 'D' -> D + 'H' -> H + 'V' -> V + _ -> N +{-# INLINE charDEG #-} + +degChar = \case + A -> 'A' + C -> 'C' + G -> 'G' + T -> 'T' + U -> 'U' + W -> 'W' + S -> 'S' + M -> 'M' + K -> 'K' + R -> 'R' + Y -> 'Y' + B -> 'B' + D -> 'D' + H -> 'H' + V -> 'V' + N -> 'N' +{-# INLINE degChar #-} + +instance Show (Letter DEG n) where + show c = [degChar c] + +degSeq :: MkPrimary p DEG n => p -> Primary DEG n +degSeq = primary + +instance MkPrimary (VU.Vector Char) DEG n where + primary = VU.map charDEG + +instance IsString [Letter DEG n] where + fromString = map charDEG + + + +-- * Conversions + +class Degenerate x where + fromDegenerate :: Char -> [x] + toDegenerate :: [x] -> Maybe Char + +instance Degenerate Char where + fromDegenerate = maybe [] id . flip lookup iupacXDNAchars + toDegenerate = flip lookup (map swap iupacXDNAchars) . nub . sort + +instance Degenerate (Letter RNA n) where + fromDegenerate 'T' = [] + fromDegenerate x = map dnaTrna $ fromDegenerate x + toDegenerate xs | xs == [R.U] = Just 'U' + | otherwise = toDegenerate $ map rnaTdna xs + +instance Degenerate (Letter DNA n) where + fromDegenerate 'U' = [] + fromDegenerate x = map charDNA $ fromDegenerate x + toDegenerate = toDegenerate . map dnaChar + +instance Degenerate (Letter XNA n) where + fromDegenerate = map charXNA . fromDegenerate + toDegenerate = toDegenerate . map xnaChar + + + +-- * Raw embeddings + +-- | list of characters, using the XNA alphabet, but degenerate chars +-- assume DNA characters. + +iupacXDNAchars :: [(Char,String)] +iupacXDNAchars = map (go . words) . lines . unpack $ iupacNucleotides where + go [[c],cs] = (c,cs) +{-# NOINLINE iupacXDNAchars #-} + +-- | Raw iupac data, embedded into the library. + +iupacNucleotides :: ByteString +iupacNucleotides = $(makeRelativeToProject "sources/iupac-nucleotides" >>= embedFile) + diff --git a/Biobase/Primary/Letter.hs b/Biobase/Primary/Letter.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Letter.hs @@ -0,0 +1,164 @@ + +-- | A newtype with an attached phantom type which allows us to encode +-- nucleotides and amino acids. Actual seqence-specific functions can be founds +-- in the appropriate modules @AA@ and @Nuc@. + +module Biobase.Primary.Letter where + +import Control.DeepSeq (NFData) +import Data.Aeson +import Data.Binary +import Data.Coerce +import Data.Data +import Data.Hashable (Hashable) +import Data.Ix (Ix(..)) +import Data.Serialize (Serialize(..)) +import Data.String (IsString(..)) +import Data.Typeable +import Data.Vector.Fusion.Stream.Monadic (map,Step(..),flatten) +import Data.Vector.Unboxed.Deriving +import GHC.Base (remInt,quotInt) +import GHC.Generics (Generic) +import Prelude hiding (map) +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as BSL +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Data.Vector.Unboxed as VU + +import Data.PrimitiveArray hiding (map) +import Data.Info + + + +-- | A 'Letter' together with its phantom type @seqTy@ encodes bio-sequences, +-- while @nameTy@ allows to specify a type-level name for a letter. + +newtype Letter (seqTy :: *) (nameTy :: k) = Letter { getLetter :: Int } + deriving (Eq,Ord,Generic,Ix,Typeable) + +-- | While @coerce@ will always work, this way restricts the change to just the +-- @nameTy@. + +changeNameTy :: Letter seqTy nameTy -> Letter seqTy newNameTy +{-# Inline changeNameTy #-} +changeNameTy = coerce + +-- | Manual @Data@ instance because @Letter@ should not show its +-- implementation. This also allows for better use of generic programming +-- downstream. + +instance (Typeable t, Typeable (Letter t n)) => Data (Letter t n) where + toConstr = mkIntegralConstr letterDataType . getLetter + gunfold _ z c = case constrRep c of + (IntConstr x) -> z (Letter $ fromIntegral x) + _ -> errorWithoutStackTrace $ "Biobase.Primary.Letter.gunfold: Constructor " + ++ show c + ++ " is not of type Letter (using Int-rep)" + dataTypeOf _ = letterDataType +letterDataType = mkDataType "Biobase.Primary.Letter" [letterConstr] +letterConstr = mkConstr letterDataType "Letter" [] Prefix + +instance Binary (Letter t n) +instance Serialize (Letter t n) + +instance NFData (Letter t n) + +type Primary t n = VU.Vector (Letter t n) + +-- | Convert 'Letter' types into character forms. @DNA@, @RNA@, and @amino +-- acid@ sequences can make use of this. Other @Letter@ types only if they +-- have single-char representations. + +class LetterChar t n where + letterChar :: Letter t n -> Char + charLetter :: Char -> Letter t n + +-- | Conversion from a large number of sequence-like inputs to primary +-- sequences. + +class MkPrimary c t n where + primary :: c -> Primary t n + +instance MkPrimary (VU.Vector Char) t n => MkPrimary String t n where + primary = primary . VU.fromList + +instance MkPrimary (VU.Vector Char) t n => MkPrimary T.Text t n where + primary = primary . VU.fromList . T.unpack + +instance MkPrimary (VU.Vector Char) t n => MkPrimary TL.Text t n where + primary = primary . VU.fromList . TL.unpack + +instance MkPrimary (VU.Vector Char) t n => MkPrimary BS.ByteString t n where + primary = primary . VU.fromList . BS.unpack + +instance MkPrimary (VU.Vector Char) t n => MkPrimary BSL.ByteString t n where + primary = primary . VU.fromList . BSL.unpack + +instance (VU.Unbox (Letter t n), IsString [Letter t n]) => IsString (VU.Vector (Letter t n)) where + fromString = VU.fromList . fromString + + + +-- *** Instances for 'Letter'. + +derivingUnbox "Letter" + [t| forall t n . Letter t n -> Int |] [| getLetter |] [| Letter |] + +instance Hashable (Letter t n) + +-- | +-- +-- TODO replace @LtLetter Int@ with more specific limits? Maybe some constants? + +instance Index (Letter l n) where + newtype LimitType (Letter l n) = LtLetter (Letter l n) + linearIndex _ (Letter i) = i + {-# Inline linearIndex #-} + fromLinearIndex _ k = Letter k + {-# Inline fromLinearIndex #-} + size (LtLetter (Letter h)) = h+1 + {-# Inline size #-} + inBounds (LtLetter h) i = zeroBound <= i && i <= h + {-# Inline inBounds #-} + zeroBound = Letter 0 + {-# Inline zeroBound #-} + zeroBound' = LtLetter zeroBound + {-# Inline zeroBound' #-} + totalSize (LtLetter (Letter k)) = [ fromIntegral k + 1 ] + {-# Inline totalSize #-} + showBound (LtLetter (Letter k)) = [ show k ] + showIndex (Letter k) = [ show k ] + +deriving instance (Bounded (Letter l n)) => Bounded (LimitType (Letter l n)) +deriving instance Eq (LimitType (Letter l n)) +deriving instance Generic (LimitType (Letter l n)) +deriving instance (Read (Letter l n)) => Read (LimitType (Letter l n)) +deriving instance (Show (Letter l n)) => Show (LimitType (Letter l n)) +deriving instance Typeable (LimitType (Letter l n)) +deriving instance Data (Letter l n) => Data (LimitType (Letter l n)) + +instance IndexStream z => IndexStream (z:.Letter l n) where + streamUp (ls:..LtLetter l) (hs:..LtLetter h) = flatten mk step $ streamUp ls hs + where mk z = return (z,l) + step (z,k) + | k > h = return $ Done + | otherwise = return $ Yield (z:.k) (z,Letter $ getLetter k +1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamUp #-} + streamDown (ls:..LtLetter l) (hs:..LtLetter h) = flatten mk step $ streamDown ls hs + where mk z = return (z,h) + step (z,k) + | k < l = return $ Done + | otherwise = return $ Yield (z:.k) (z,Letter $ getLetter k -1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamDown #-} + +instance IndexStream (Letter l n) where + streamUp l h = map (\(Z:.k) -> k) $ streamUp (ZZ:..l) (ZZ:..h) + streamDown l h = map (\(Z:.k) -> k) $ streamDown (ZZ:..l) (ZZ:..h) + {-# Inline streamUp #-} + {-# Inline streamDown #-} + diff --git a/Biobase/Primary/Nuc.hs b/Biobase/Primary/Nuc.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Nuc.hs @@ -0,0 +1,28 @@ + +-- | The primary structure: interface to efficient encoding of RNA and DNA +-- sequences. The design aims toward the 'vector' library and repa. In +-- particular, everything is strict; if you want to stream full genomes, use +-- 'text' or lazy 'bytestring's instead and cast to Biobase.Primary definitions +-- only at the last moment. +-- +-- Degenerate encoding can be found in the @IUPAC@ module. +-- +-- TODO enable OverloadedLists + +module Biobase.Primary.Nuc + ( module Biobase.Primary.Letter + , module Biobase.Primary.Nuc.Conversion + , module Biobase.Primary.Nuc.DNA + , module Biobase.Primary.Nuc.RNA + , module Biobase.Primary.Nuc.XNA + ) where + +import Biobase.Primary.Letter +import Biobase.Primary.Nuc.Conversion +import Biobase.Primary.Nuc.DNA hiding (A,C,G,T,N) +import Biobase.Primary.Nuc.RNA hiding (A,C,G,U,N) +import Biobase.Primary.Nuc.XNA hiding (A,C,G,T,U,N) +import qualified Biobase.Primary.Nuc.DNA as D +import qualified Biobase.Primary.Nuc.RNA as R +import qualified Biobase.Primary.Nuc.XNA as X + diff --git a/Biobase/Primary/Nuc/Conversion.hs b/Biobase/Primary/Nuc/Conversion.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Nuc/Conversion.hs @@ -0,0 +1,179 @@ + +-- {-# Language CPP #-} +-- +-- #if __GLASGOW_HASKELL__ < 710 +-- {-# LANGUAGE OverlappingInstances #-} +-- #endif + +-- | Convert between different nucleotide representations + +module Biobase.Primary.Nuc.Conversion where + +import Control.Lens (iso, from) +import qualified Data.Vector.Unboxed as VU + +import Biobase.Types.BioSequence (Transcribe(..),RNA,DNA) + +import Biobase.Primary.Letter (Letter(..), Primary) +import qualified Biobase.Primary.Nuc.DNA as D +import qualified Biobase.Primary.Nuc.RNA as R +import qualified Biobase.Primary.Nuc.XNA as X + + + +-- * Single-character translations. + +-- | Transform RNA to DNA. That means change @U@ to @T@ and keep the other +-- characters as is. + +rnaTdna = \case + R.A -> D.A + R.C -> D.C + R.G -> D.G + R.U -> D.T + _ -> D.N +{-# INLINE rnaTdna #-} + +-- | Transform DNA to RNA. That means change @T@ to @U@ and keep the other +-- characters as is. + +dnaTrna = \case + D.A -> R.A + D.C -> R.C + D.G -> R.G + D.T -> R.U + _ -> R.N +{-# INLINE dnaTrna #-} + +-- | Generalize an RNA character to a XNA character. + +rnaGxna = \case + R.A -> X.A + R.C -> X.C + R.G -> X.G + R.U -> X.U + _ -> X.N +{-# INLINE rnaGxna #-} + +-- | Generalize a DNA character to a XNA character. + +dnaGxna = \case + D.A -> X.A + D.C -> X.C + D.G -> X.G + D.T -> X.T + _ -> X.N +{-# INLINE dnaGxna #-} + +-- | Specialize XNA to RNA, @T@ becomes @N@. + +xnaSrna = \case + X.A -> R.A + X.C -> R.C + X.G -> R.G + X.U -> R.U + _ -> R.N +{-# INLINE xnaSrna #-} + +-- | Specialize XNA to DNA, @U@ becomes @N@. + +xnaSdna = \case + X.A -> D.A + X.C -> D.C + X.G -> D.G + X.T -> D.T + _ -> D.N +{-# INLINE xnaSdna #-} + + + +-- ** Transcription between RNA and DNA. Both on the individual sequence level, +-- and on the level of primary sequence data. + +instance Transcribe (Letter RNA n) where + type TranscribeTo (Letter RNA n) = Letter DNA n + transcribe = iso rnaTdna dnaTrna + {-# Inline transcribe #-} + +instance Transcribe (Letter DNA n) where + type TranscribeTo (Letter DNA n) = Letter RNA n + transcribe = from transcribe + {-# Inline transcribe #-} + +instance Transcribe (Primary RNA n) where + type TranscribeTo (Primary RNA n) = Primary DNA n + transcribe = iso (VU.map rnaTdna) (VU.map dnaTrna) + {-# Inline transcribe #-} + +instance Transcribe (Primary DNA n) where + type TranscribeTo (Primary DNA n) = Primary RNA n + transcribe = iso (VU.map dnaTrna) (VU.map rnaTdna) + {-# Inline transcribe #-} + + +-- TODO to be removed soon + +---- * Reverse-complement of characters. +-- +---- | Produce the complement of a RNA or DNA sequence. Does intentionally +---- not work for XNA sequences as it is not possible to uniquely translate +---- @A@ into either @U@ or @T@. +-- +--class Complement s t where +-- complement :: s -> t +-- +---- | To 'transcribe' a DNA sequence into RNA we reverse the complement of +---- the sequence. +-- +--transcribe :: Primary D.DNA -> Primary R.RNA +--transcribe = VU.reverse . complement +-- +--instance Complement (Letter R.RNA) (Letter R.RNA) where +-- complement = \case +-- R.A -> R.U +-- R.C -> R.G +-- R.G -> R.C +-- R.U -> R.A +-- R.N -> R.N +-- +--instance Complement (Letter D.DNA) (Letter D.DNA) where +-- complement = \case +-- D.A -> D.T +-- D.C -> D.G +-- D.G -> D.C +-- D.T -> D.A +-- D.N -> D.N +-- +--instance Complement (Letter D.DNA) (Letter R.RNA) where +-- complement = \case +-- D.A -> R.U +-- D.C -> R.G +-- D.G -> R.C +-- D.T -> R.A +-- D.N -> R.N +-- +--instance Complement (Letter R.RNA) (Letter D.DNA) where +-- complement = \case +-- R.A -> D.T +-- R.C -> D.G +-- R.G -> D.C +-- R.U -> D.A +-- R.N -> D.N +-- +-- #if __GLASGOW_HASKELL__ >= 710 +-- instance {-# OVERLAPPING #-} +-- #else +-- instance +-- #endif +-- ( Complement s t, VU.Unbox s, VU.Unbox t) +-- => Complement (VU.Vector s) (VU.Vector t) +-- where complement = VU.map complement +-- +-- #if __GLASGOW_HASKELL__ >= 710 +-- instance {-# Overlappable #-} +-- #else +-- instance +-- #endif +-- ( Complement s t, Functor f) => Complement (f s) (f t) +-- where complement = fmap complement + diff --git a/Biobase/Primary/Nuc/DNA.hs b/Biobase/Primary/Nuc/DNA.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Nuc/DNA.hs @@ -0,0 +1,101 @@ + +module Biobase.Primary.Nuc.DNA where + +import Control.Category ((>>>)) +import Control.Lens (Iso', iso) +import Data.Aeson +import Data.Char (toUpper) +import Data.Ix (Ix(..)) +import Data.Primitive.Types +import Data.String +import Data.Tuple (swap) +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as BSL +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Unboxed as VU + +import Biobase.Types.BioSequence (DNA) + +import Biobase.Primary.Bounds +import Biobase.Primary.Letter + + + +-- Single-character names for nucleotides. + +pattern A = Letter 0 :: Letter DNA n +pattern C = Letter 1 :: Letter DNA n +pattern G = Letter 2 :: Letter DNA n +pattern T = Letter 3 :: Letter DNA n +pattern N = Letter 4 :: Letter DNA n + +instance Enum (Letter DNA n) where + succ N = error "succ/N:DNA" + succ (Letter x) = Letter $ x+1 + pred A = error "pred/A:DNA" + pred (Letter x) = Letter $ x-1 + toEnum k | k>=0 && k<=4 = Letter k + toEnum k = error $ "toEnum/Letter DNA " ++ show k + fromEnum (Letter k) = k + +instance LetterChar DNA n where + letterChar = dnaChar + charLetter = charDNA + +--instance (LetterChar DNA) => ToJSON (Primary DNA) where +-- toJSON = toJSON . VU.toList . VU.map letterChar +-- +--instance (MkPrimary (VU.Vector Char) DNA) => FromJSON (Primary DNA) where +-- parseJSON = fmap (primary :: String -> Primary DNA) . parseJSON + +acgt :: [Letter DNA n] +acgt = [A .. T] + +charDNA = toUpper >>> \case + 'A' -> A + 'C' -> C + 'G' -> G + 'T' -> T + _ -> N +{-# INLINE charDNA #-} + +dnaChar = \case + A -> 'A' + C -> 'C' + G -> 'G' + T -> 'T' + N -> 'N' +{-# INLINE dnaChar #-} + +-- | An isomorphism from 'Char' to 'Letter DNA'. This assumes that the +-- underlying @Char@s actually represent a DNA sequence. This allows typesafe +-- modification of DNA sequences since only @[A,C,G,T,N]@ are allowed. + +cdna ∷ Iso' Char (Letter DNA n) +cdna = iso charDNA dnaChar + +instance Show (Letter DNA n) where + show c = [dnaChar c] + +instance Read (Letter DNA n) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | otherwise = [(charDNA x, xs)] + +dnaSeq :: MkPrimary p DNA n => p -> Primary DNA n +dnaSeq = primary + +instance Bounded (Letter DNA n) where + minBound = A + maxBound = N + +instance MkPrimary (VU.Vector Char) DNA n where + primary = VU.map charDNA + +instance IsString [Letter DNA n] where + fromString = map charDNA + diff --git a/Biobase/Primary/Nuc/RNA.hs b/Biobase/Primary/Nuc/RNA.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Nuc/RNA.hs @@ -0,0 +1,118 @@ + +module Biobase.Primary.Nuc.RNA where + +import Control.Category ((>>>)) +import Control.Lens (Iso', iso) +import Data.Aeson +import Data.Char (toUpper) +import Data.Data +import Data.Ix (Ix(..)) +import Data.Primitive.Types +import Data.String +import Data.Tuple (swap) +import Data.Typeable +import qualified Data.ByteString.Builder as BB +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as BSL +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Unboxed as VU + +import Biobase.Types.BioSequence (RNA) + +import Biobase.Primary.Bounds +import Biobase.Primary.Letter + + + +pattern A = Letter 0 ∷ Letter RNA n +pattern C = Letter 1 ∷ Letter RNA n +pattern G = Letter 2 ∷ Letter RNA n +pattern U = Letter 3 ∷ Letter RNA n +pattern N = Letter 4 ∷ Letter RNA n + +instance Bounded (Letter RNA n) where + minBound = A + maxBound = N + +instance Enum (Letter RNA n) where + succ N = error "succ/N:RNA" + succ (Letter x) = Letter $ x+1 + pred A = error "pred/A:RNA" + pred (Letter x) = Letter $ x-1 + toEnum k | k>=0 && k<=4 = Letter k + toEnum k = error $ "toEnum/Letter RNA " ++ show k + fromEnum (Letter k) = k + +instance LetterChar RNA n where + letterChar = rnaChar + charLetter = charRNA + +instance ToJSON (Letter RNA n) where + toJSON = toJSON . letterChar + +instance FromJSON (Letter RNA n) where + parseJSON = fmap charLetter . parseJSON + +-- We encode 'Primary RNA' directly as a string. +-- +-- TODO we can't anymore, because this is not a newtype, just a type. + +--instance ToJSON (Primary RNA) where +-- toJSON = toJSON . VU.toList . VU.map letterChar +-- +--instance FromJSON (Primary RNA) where +-- parseJSON = fmap (primary ∷ String → Primary RNA) . parseJSON + + +acgu ∷ [Letter RNA n] +acgu = [A .. U] + +charRNA = toUpper >>> \case + 'A' -> A + 'C' -> C + 'G' -> G + 'U' -> U + _ -> N +{-# INLINE charRNA #-} + +rnaChar = \case + A -> 'A' + C -> 'C' + G -> 'G' + U -> 'U' + N -> 'N' + _ -> '\9888' +{-# INLINE rnaChar #-} + +-- | An isomorphism from 'Char' to 'Letter RNA'. This assumes that the +-- underlying @Char@s actually represent an RNA sequence. This allows typesafe +-- modification of RNA sequences since only @[A,C,G,U,N]@ are allowed. + +crna ∷ Iso' Char (Letter RNA n) +crna = iso charRNA rnaChar + +instance Show (Letter RNA n) where + show c = [rnaChar c] + +instance Read (Letter RNA n) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | otherwise = [(charRNA x, xs)] + +rnaSeq ∷ MkPrimary p RNA n ⇒ p → Primary RNA n +rnaSeq = primary + +instance MkPrimary (VU.Vector Char) RNA n where + primary = VU.map charRNA + +instance IsString [Letter RNA n] where + fromString = map charRNA + +viennaPairs = [ (C,G), (G,C), (G,U), (U,G), (A,U), (U,A) ] +viennaPairsNN = viennaPairs ++ [ (N,N) ] + + diff --git a/Biobase/Primary/Nuc/XNA.hs b/Biobase/Primary/Nuc/XNA.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Nuc/XNA.hs @@ -0,0 +1,93 @@ + +module Biobase.Primary.Nuc.XNA where + +import Data.Aeson +import Data.Char (toUpper) +import Data.Ix (Ix(..)) +import Data.Primitive.Types +import Data.String +import Data.Tuple (swap) +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as BSL +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Unboxed as VU +import Control.Category ((>>>)) + +import Biobase.Types.BioSequence + +import Biobase.Primary.Bounds +import Biobase.Primary.Letter + + + +-- | Combine both, RNA and DNA. + +pattern A = Letter 0 :: Letter XNA n +pattern C = Letter 1 :: Letter XNA n +pattern G = Letter 2 :: Letter XNA n +pattern T = Letter 3 :: Letter XNA n +pattern U = Letter 4 :: Letter XNA n +pattern N = Letter 5 :: Letter XNA n + +instance Bounded (Letter XNA n) where + minBound = A + maxBound = N + +instance Enum (Letter XNA n) where + succ N = error "succ/N:XNA" + succ (Letter x) = Letter $ x+1 + pred A = error "pred/A:XNA" + pred (Letter x) = Letter $ x-1 + toEnum k | k>=0 && k<=5 = Letter k + toEnum k = error $ "toEnum/Letter XNA " ++ show k + fromEnum (Letter k) = k + +instance LetterChar XNA n where + letterChar = xnaChar + charLetter = charXNA + +--instance (LetterChar XNA) => ToJSON (Primary XNA) where +-- toJSON = toJSON . VU.toList . VU.map letterChar +-- +--instance (MkPrimary (VU.Vector Char) XNA) => FromJSON (Primary XNA) where +-- parseJSON = fmap (primary :: String -> Primary XNA) . parseJSON + +charXNA = toUpper >>> \case + 'A' -> A + 'C' -> C + 'G' -> G + 'T' -> T + 'U' -> U + _ -> N +{-# INLINE charXNA #-} + +xnaChar = \case + A -> 'A' + C -> 'C' + G -> 'G' + T -> 'T' + U -> 'U' + N -> 'N' +{-# INLINE xnaChar #-} + +instance Show (Letter XNA n) where + show c = [xnaChar c] + +instance Read (Letter XNA n) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | otherwise = [(charXNA x, xs)] + +xnaSeq :: MkPrimary p XNA n => p -> Primary XNA n +xnaSeq = primary + +instance MkPrimary (VU.Vector Char) XNA n where + primary = VU.map charXNA + +instance IsString [Letter XNA n] where + fromString = map charXNA + diff --git a/Biobase/Primary/Pretty.hs b/Biobase/Primary/Pretty.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Pretty.hs @@ -0,0 +1,30 @@ + +-- | Wrapper newtype to simplify pretty and short encoding of primary +-- sequences. + +module Biobase.Primary.Pretty where + +import Data.Aeson +import qualified Data.Vector.Unboxed as VU +import qualified Data.Vector as V +import qualified Data.Vector.Storable as VS +import qualified Data.Text as T + +import Biobase.Primary.Letter + + + +newtype Pretty f a = Pretty { getPretty :: f a } + +instance (LetterChar x n) => ToJSON (Pretty VU.Vector (Letter x n)) where + toJSON = String . T.pack . map letterChar . VU.toList . getPretty + +instance (LetterChar x n) => ToJSON (Pretty V.Vector (Letter x n)) where + toJSON = String . T.pack . map letterChar . V.toList . getPretty + +instance (LetterChar x n, VS.Storable (Letter x n)) => ToJSON (Pretty VS.Vector (Letter x n)) where + toJSON = String . T.pack . map letterChar . VS.toList . getPretty + +instance (LetterChar x n) => ToJSON (Pretty [] (Letter x n)) where + toJSON = String . T.pack . map letterChar . getPretty + diff --git a/Biobase/Primary/Trans.hs b/Biobase/Primary/Trans.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Trans.hs @@ -0,0 +1,79 @@ + +-- | This module provides functionality for translation between nucleotides +-- and amino acids. +-- +-- NOTE 'aaDNAseq' is lossy. Might be a good idea to consider something +-- more involved? +-- +-- TODO we need different functions, depending on if we have a part of +-- a genome in @DNA@ form, or some messenger @RNA@. It'll probably also be +-- useful to return @Either@, with @Left@ indicating error like partially +-- translated sequence due to intermediate stop codons, or so. +-- +-- TODO 'dnaAAseq' and 'aaDNAseq' can be nicely optimized using 'flatten' +-- and friends. + +module Biobase.Primary.Trans where + +import Control.Lens +import Control.Arrow ((***)) +import Data.ByteString.Char8 (ByteString,unpack) +import Data.Map.Strict (Map) +import Data.Tuple (swap) +import qualified Data.Map.Strict as M +import qualified Data.Vector.Unboxed as VU + +import Biobase.Types.BioSequence +import Biobase.Types.Codon + +import Biobase.Primary.AA +import Biobase.Primary.Nuc +import Biobase.Primary.Letter +import Biobase.GeneticCodes.Translation +import Biobase.GeneticCodes.Types + + + +-- | Transform translation tables into the @Letter DNA/Letter AA@ format. + +letterTranslationTable :: TranslationTable Char Char -> TranslationTable (Letter DNA n) (Letter AA n) +letterTranslationTable tbl = TranslationTable + { _codonToAminoAcid = M.fromList . map (ftriplet *** felement) . M.toList $ tbl^.codonToAminoAcid + , _aminoAcidtoCodons = M.fromList . map (charAA *** map felement) . M.toList $ tbl^.aminoAcidtoCodons + , _tableID = tbl^.tableID + , _tableName = tbl^.tableName + } where ftriplet :: Codon Char -> Codon (Letter DNA n) + ftriplet = over each charDNA + felement :: TranslationElement Char Char -> TranslationElement (Letter DNA n) (Letter AA n) + felement = over (baseCodon.each) charDNA . over aminoAcid charAA + +instance Translation (Codon (Letter DNA n)) where + type TargetType (Codon (Letter DNA n)) = Letter AA n + type CodonType (Codon (Letter DNA n)) = Letter DNA n + type AAType (Codon (Letter DNA n)) = Letter AA n + translate tbl t = maybe Unknown _aminoAcid $ M.lookup t (tbl^.codonToAminoAcid) + {-# Inline translate #-} + translateAllFrames = translate + {-# Inline translateAllFrames #-} + +instance Translation (Primary DNA n) where + type TargetType (Primary DNA n) = Primary AA n + type CodonType (Primary DNA n) = Letter DNA n + type AAType (Primary DNA n) = Letter AA n + -- | + -- + -- TODO we could consider returning @Nothing@ in case the input is not + -- power-of-three. + translate tbl xs = VU.unfoldrN (VU.length xs `div` 3) go xs + where go (VU.splitAt 3 -> (hs,ts)) + | VU.length hs < 3 = Nothing + | otherwise = Just (aa,ts) + where [a,b,c] = VU.toList hs + aa = translate tbl $ Codon a b c + {-# Inline translate #-} + translateAllFrames tbl xs = VU.unfoldrN (VU.length xs) go 0 + where go 0 = Just (Undef,1) + go 1 = Just (Undef,2) + go k = Just (translate tbl $ Codon (xs VU.! (k-2)) (xs VU.! (k-1)) (xs VU.! k), k+1) + {-# Inlinable translateAllFrames #-} + diff --git a/Biobase/Primary/Unknown.hs b/Biobase/Primary/Unknown.hs new file mode 100644 --- /dev/null +++ b/Biobase/Primary/Unknown.hs @@ -0,0 +1,82 @@ + +-- | A 'Letter' with unknown annotation. We sometimes want to encode that we +-- are dealing with @Letter@s in an alphabet, but we do not want to commit to a +-- certain alphabet (just yet). +-- +-- This module allows us to make explicit that we do not know the specific +-- alphabet type yet. +-- +-- One should NEVER blindly coerce, since the order and limits of @Letter@'s +-- might well be different. + +module Biobase.Primary.Unknown where + +import Data.Aeson +import Control.Applicative ((<$>)) +import Control.Arrow ((***),first) +import Data.Hashable +import Data.Ix (Ix(..)) +import Data.Map.Strict (Map) +import Data.Primitive.Types +import Data.Tuple (swap) +import Data.Vector.Unboxed.Deriving +import Debug.Trace +import GHC.Base (remInt,quotInt) +import GHC.Generics (Generic) +import GHC.Read +import qualified Data.Bijection.Map as B +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as BSL +import qualified Data.Map.Strict as M +import qualified Data.Text as T +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Unboxed as VU +import qualified Text.ParserCombinators.ReadPrec as RP +import qualified Text.Read.Lex as Lex + +import Biobase.Primary.Letter + + + +-- | @Unknown@ phantom type. + +data Unknown + + + +-- | Creating an unknown letter. + +unk ∷ Int → Letter Unknown n +unk = Letter + + + +-- *** instances + +instance Show (Letter Unknown n) where + show (Letter i) = "U " ++ show i + +instance Read (Letter Unknown n) where + readPrec = parens $ do + Lex.Ident u <- lexP + case u of + "U" → unk <$> readPrec + _ → RP.pfail + +instance Enum (Letter Unknown n) where + succ (Letter x) = Letter $ x+1 + pred (Letter x) = Letter $ x-1 + toEnum = Letter + fromEnum = getLetter + +instance MkPrimary (VU.Vector Int) Unknown n where + primary = VU.map Letter + {-# Inline primary #-} + +instance ToJSON (Letter Unknown n) where + toJSON = toJSON . getLetter + +instance FromJSON (Letter Unknown n) where + parseJSON = fmap Letter . parseJSON + diff --git a/Biobase/Secondary.hs b/Biobase/Secondary.hs --- a/Biobase/Secondary.hs +++ b/Biobase/Secondary.hs @@ -1,348 +1,19 @@ -{-# LANGUAGE PackageImports #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE OverlappingInstances #-} -{-# LANGUAGE StandaloneDeriving #-} --- | Secondary structure: define basepairs as Int-tuples, the three edges, a --- nucleotide can use for pairing and the cis/trans isomerism. Both edges and --- cis/trans come with a tag for "unknown". --- --- TODO set ext-annotations to be (isomerism,edge,edge) and have a asString --- instance to read "cWW" "tSH" and other notation. - -module Biobase.Secondary where - -import "PrimitiveArray" Data.Array.Repa.Index -import "PrimitiveArray" Data.Array.Repa.Shape -import Data.Char (toLower, toUpper) -import Data.Ix (Ix(..)) -import Data.List as L -import Data.Primitive.Types -import Data.Tuple (swap) -import GHC.Base (remInt,quotInt) -import qualified Data.Vector.Generic as VG -import qualified Data.Vector.Generic.Mutable as VGM -import qualified Data.Vector.Unboxed as VU - -import Biobase.Primary -import Biobase.Primary.Bounds - - - --- | Easy reading of a three-Char string into a triple. - -threeChar :: String -> ExtPairAnnotation -threeChar s@[c,x,y] - | Just c' <- L.lookup (toLower c) charCTList - , Just x' <- L.lookup (toUpper x) charEdgeList - , Just y' <- L.lookup (toUpper y) charEdgeList - = (c',x',y') - | map toLower s == "bif" = (unknownCT,unknownEdge,unknownEdge) - | otherwise = error $ "can't convert string: " ++ s - --- | Each nucleotide in a pair may be paired using one of three edges: --- watson-crik, sugar, or hoogsteen. - -newtype Edge = Edge {unEdge :: Int} - deriving (Eq,Ord,Ix) - -instance (Shape sh,Show sh) => Shape (sh :. Edge) where - rank (sh:._) = rank sh + 1 - zeroDim = zeroDim:.Edge 0 - unitDim = unitDim:.Edge 1 -- TODO does this one make sense? - intersectDim (sh1:.n1) (sh2:.n2) = intersectDim sh1 sh2 :. min n1 n2 - addDim (sh1:.Edge n1) (sh2:.Edge n2) = addDim sh1 sh2 :. Edge (n1+n2) -- TODO will not necessarily yield a valid Edge - size (sh1:.Edge n) = size sh1 * n - sizeIsValid (sh1:.Edge n) = sizeIsValid (sh1:.n) - toIndex (sh1:.Edge sh2) (sh1':.Edge sh2') = toIndex (sh1:.sh2) (sh1':.sh2') - fromIndex (ds:.Edge d) n = fromIndex ds (n `quotInt` d) :. Edge r where - r | rank ds == 0 = n - | otherwise = n `remInt` d - inShapeRange (sh1:.n1) (sh2:.n2) (idx:.i) = i>=n1 && i error "empty list in shapeOfList/Primary" - x:xs -> shapeOfList xs :. Edge x - deepSeq (sh:.n) x = deepSeq sh (n `seq` x) - {-# INLINE rank #-} - {-# INLINE zeroDim #-} - {-# INLINE unitDim #-} - {-# INLINE intersectDim #-} - {-# INLINE addDim #-} - {-# INLINE size #-} - {-# INLINE sizeIsValid #-} - {-# INLINE toIndex #-} - {-# INLINE fromIndex #-} - {-# INLINE inShapeRange #-} - {-# INLINE listOfShape #-} - {-# INLINE shapeOfList #-} - {-# INLINE deepSeq #-} - -(wc : sugar : hoogsteen : unknownEdge : edgeUndefined : _) = map Edge [0..] - -charEdgeList = - [ ('W',wc) - , ('S',sugar) - , ('H',hoogsteen) - , ('?',unknownEdge) - ] - -edgeCharList = map swap charEdgeList - --- | Human-readable Show instance. - -instance Show Edge where - show k - | Just v <- k `lookup` edgeCharList = [v] - | otherwise = "?" - --- | Human-readable Read instance. - -instance Read Edge where - readsPrec p [] = [] - readsPrec p (x:xs) - | x ==' ' = readsPrec p xs - | Just n <- x `lookup` charEdgeList = [(n,xs)] - | otherwise = [] - --- | Nucleotides in a pairing may be in the cis(==?) or trans(==?) state. - -newtype CTisomerism = CT {unCT :: Int} - deriving (Eq,Ord,Ix) - -instance (Shape sh,Show sh) => Shape (sh :. CTisomerism) where - rank (sh:._) = rank sh + 1 - zeroDim = zeroDim:.CT 0 - unitDim = unitDim:.CT 1 -- TODO does this one make sense? - intersectDim (sh1:.n1) (sh2:.n2) = intersectDim sh1 sh2 :. min n1 n2 - addDim (sh1:.CT n1) (sh2:.CT n2) = addDim sh1 sh2 :. CT (n1+n2) -- TODO will not necessarily yield a valid CT - size (sh1:.CT n) = size sh1 * n - sizeIsValid (sh1:.CT n) = sizeIsValid (sh1:.n) - toIndex (sh1:.CT sh2) (sh1':.CT sh2') = toIndex (sh1:.sh2) (sh1':.sh2') - fromIndex (ds:.CT d) n = fromIndex ds (n `quotInt` d) :. CT r where - r | rank ds == 0 = n - | otherwise = n `remInt` d - inShapeRange (sh1:.n1) (sh2:.n2) (idx:.i) = i>=n1 && i error "empty list in shapeOfList/Primary" - x:xs -> shapeOfList xs :. CT x - deepSeq (sh:.n) x = deepSeq sh (n `seq` x) - {-# INLINE rank #-} - {-# INLINE zeroDim #-} - {-# INLINE unitDim #-} - {-# INLINE intersectDim #-} - {-# INLINE addDim #-} - {-# INLINE size #-} - {-# INLINE sizeIsValid #-} - {-# INLINE toIndex #-} - {-# INLINE fromIndex #-} - {-# INLINE inShapeRange #-} - {-# INLINE listOfShape #-} - {-# INLINE shapeOfList #-} - {-# INLINE deepSeq #-} - -(cis : trans : unknownCT : undefinedCT : _) = map CT [0..] -antiCT = undefined -paraCT = undefined - -charCTList = - [ ('c',cis) - , ('t',trans) - , ('?',unknownCT) - -- TODO antiCT, paraCT - -- TODO '?' type (??? could denote bifurcation) - ] - -ctCharList = map swap charCTList - --- | Human-readable Show instance. - -instance Show CTisomerism where - show k - | Just v <- k `lookup` ctCharList = [v] - | otherwise = "?" - --- | Human-readable Read instance. - -instance Read CTisomerism where - readsPrec p [] = [] - readsPrec p (x:xs) - | x ==' ' = readsPrec p xs - | Just n <- x `lookup` charCTList = [(n,xs)] - | otherwise = [] - - - --- * Instances - --- ** Instances for 'Edge' - -deriving instance Prim Edge -deriving instance VGM.MVector VU.MVector Edge -deriving instance VG.Vector VU.Vector Edge -deriving instance VU.Unbox Edge - -instance Bounded Edge where - minBound = wc - maxBound = unknownEdge - -instance Bounds Edge where - minNormal = wc - maxNormal = wc - minExtended = wc - maxExtended = hoogsteen - -instance Enum Edge where - toEnum = Edge - fromEnum = unEdge - --- ** Instances for 'CTisomerism' - -deriving instance Prim CTisomerism -deriving instance VGM.MVector VU.MVector CTisomerism -deriving instance VG.Vector VU.Vector CTisomerism -deriving instance VU.Unbox CTisomerism - -instance Bounded CTisomerism where - minBound = cis - maxBound = unknownCT - -instance Bounds CTisomerism where - minNormal = cis - maxNormal = cis - minExtended = cis - maxExtended = trans - -instance Enum CTisomerism where - toEnum = CT - fromEnum = unCT - - - --- * Types - --- | A basepair is simply a pair of Ints which are 0-indexing a sequence. --- --- TODO storable vector, newtype, peek/poke? - -type PairIdx = (Int,Int) - --- | A pair as a tuple containing 'Nuc's. - -type Pair = (Nuc,Nuc) - --- | Annotation for a basepair. - -type ExtPairAnnotation = (CTisomerism,Edge,Edge) - --- | An extended basepair is a basepair, annotated with edge and CTisomerism. - -type ExtPairIdx = (PairIdx,ExtPairAnnotation) - --- | An extended basepair, with nucleotides an annotation. - -type ExtPair = (Pair,ExtPairAnnotation) - - - --- * little helpers - -cWW = (cis,wc,wc) -cWS = (cis,wc,sugar) -cWH = (cis,wc,hoogsteen) -cSW = (cis,sugar,wc) -cSS = (cis,sugar,sugar) -cSH = (cis,sugar,hoogsteen) -cHW = (cis,hoogsteen,wc) -cHS = (cis,hoogsteen,sugar) -cHH = (cis,hoogsteen,hoogsteen) -tWW = (trans,wc,wc) -tWS = (trans,wc,sugar) -tWH = (trans,wc,hoogsteen) -tSW = (trans,sugar,wc) -tSS = (trans,sugar,sugar) -tSH = (trans,sugar,hoogsteen) -tHW = (trans,hoogsteen,wc) -tHS = (trans,hoogsteen,sugar) -tHH = (trans,hoogsteen,hoogsteen) - - --- * special show instances - --- | This one requires ghc head --- --- TODO maybe newtype this triple? - ---instance Show (CTisomerism,Edge,Edge) where --- show (ct,eI,eJ) = concat [show ct, show eI, show eJ] - - - --- * tuple-like selection - --- | Selection of nucleotides and/or type classes independent of which type we --- are looking at. - -class BaseSelect a b | a -> b where - -- | select first index or nucleotide - baseL :: a -> b - -- | select second index or nucleotide - baseR :: a -> b - -- | select both nucleotides as pair - baseP :: a -> (b,b) - -- | select basepair type if existing or return default cWW - baseT :: a -> ExtPairAnnotation - -- | update first index or nucleotide - updL :: b -> a -> a - -- | update second index or nucleotide - updR :: b -> a -> a - -- | update complete pair - updP :: (b,b) -> a -> a - -- | update basepair type, error if not possible due to type a - updT :: ExtPairAnnotation -> a -> a - --- | extended pairtype annotation given - -instance BaseSelect ((a,a),ExtPairAnnotation) a where - baseL ((a,_),_) = a - baseR ((_,b),_) = b - baseP (lr ,_) = lr - baseT (_,t) = t - updL n ((_,y),t) = ((n,y),t) - updR n ((x,_),t) = ((x,n),t) - updP n (_,t) = (n,t) - updT n (xy,_) = (xy,n) - {-# INLINE baseL #-} - {-# INLINE baseR #-} - {-# INLINE baseP #-} - {-# INLINE baseT #-} - {-# INLINE updL #-} - {-# INLINE updR #-} - {-# INLINE updP #-} - {-# INLINE updT #-} +module Biobase.Secondary +-- ( module Biobase.Secondary.Basepair +-- ( module Biobase.Secondary.Constraint + ( module Biobase.Secondary.Diagrams + , module Biobase.Secondary.Isostericity + , module Biobase.Secondary.Pseudoknots + , module Biobase.Secondary.Structure +-- , module Biobase.Secondary.Vienna + ) where --- | simple cis/wc-wc basepairs +--import Biobase.Secondary.Basepair +--import Biobase.Secondary.Constraint +import Biobase.Secondary.Diagrams +import Biobase.Secondary.Isostericity +import Biobase.Secondary.Pseudoknots +import Biobase.Secondary.Structure +--import Biobase.Secondary.Vienna -instance BaseSelect (a,a) a where - baseL (a,_) = a - baseR (_,a) = a - baseP = id - baseT _ = cWW - updL n (_,y) = (n,y) - updR n (x,_) = (x,n) - updP n _ = n - updT n xy = if n==cWW then xy else error $ "updT on standard pairs can not update to: " ++ show n - {-# INLINE baseL #-} - {-# INLINE baseR #-} - {-# INLINE baseP #-} - {-# INLINE baseT #-} - {-# INLINE updL #-} - {-# INLINE updR #-} - {-# INLINE updP #-} - {-# INLINE updT #-} diff --git a/Biobase/Secondary/Basepair.hs b/Biobase/Secondary/Basepair.hs new file mode 100644 --- /dev/null +++ b/Biobase/Secondary/Basepair.hs @@ -0,0 +1,260 @@ + +-- | Secondary structure: define basepairs as Int-tuples, the three edges, a +-- nucleotide can use for pairing and the cis/trans isomerism. Both edges and +-- cis/trans come with a tag for "unknown". +-- +-- Since we often want to make "pairedness" explicit, we have a newtype for +-- this as well. +-- +-- TODO set ext-annotations to be (isomerism,edge,edge) and have a asString +-- instance to read "cWW" "tSH" and other notation. + +module Biobase.Secondary.Basepair where + +import Data.Aeson +import Data.Binary +import Data.Char (toLower, toUpper) +import Data.Ix (Ix(..)) +import Data.List as L +import Data.Primitive.Types +import Data.Serialize (Serialize) +import Data.Tuple (swap) +import Data.Vector.Fusion.Stream.Monadic (map,Step(..),flatten) +import Data.Vector.Unboxed.Deriving +import GHC.Base (remInt,quotInt) +import GHC.Generics +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Unboxed as VU +import Text.Read + +import Biobase.Types.BioSequence +import Data.PrimitiveArray hiding (Complement(..),map) + +import Biobase.Primary +import Biobase.Primary.Nuc.RNA +import Biobase.Primary.Nuc + + + +-- * Newtype for efficient basepair encoding. + +-- | Encode a base pair as a single @Int@. + +newtype Basepair = BP { getBP :: Int } + deriving (Eq,Ord,Ix,Generic) + +derivingUnbox "Basepair" + [t| Basepair -> Int |] [| getBP |] [| BP |] + +instance Binary Basepair +instance Serialize Basepair +instance FromJSON Basepair +instance ToJSON Basepair + +instance Index Basepair where + newtype LimitType Basepair = LtBP Basepair + +instance IndexStream z => IndexStream (z:.Basepair) where + streamUp (ls:..LtBP (BP l)) (hs:..LtBP (BP h)) = flatten mk step $ streamUp ls hs + where mk z = return (z,l) + step (z,k) + | k > h = return $ Done + | otherwise = return $ Yield (z:.BP k) (z,k+1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamUp #-} + streamDown (ls:..LtBP (BP l)) (hs:..LtBP (BP h)) = flatten mk step $ streamDown ls hs + where mk z = return (z,h) + step (z,k) + | k < l = return $ Done + | otherwise = return $ Yield (z:.BP k) (z,k-1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamDown #-} + +instance IndexStream Basepair + +pattern AA = BP 0 +pattern AC = BP 1 +pattern AG = BP 2 +pattern AU = BP 3 +pattern CA = BP 4 +pattern CC = BP 5 +pattern CG = BP 6 +pattern CU = BP 7 +pattern GA = BP 8 +pattern GC = BP 9 +pattern GG = BP 10 +pattern GU = BP 11 +pattern UA = BP 12 +pattern UC = BP 13 +pattern UG = BP 14 +pattern UU = BP 15 +pattern NS = BP 16 +pattern NoBP = BP 17 + +{- +class MkBasepair a where + mkBasepair :: a -> Basepair + fromBasepair :: Basepair -> a + +-- | If we get a "legal" base pair, we just create it, all other +-- combinations yield 'NoBP'. Non-standard base pairs have to be created +-- explicitly using @NS@. When going back to @a@, non-standard and no pair +-- yield @(N,N)@. + +instance MkBasepair (Letter RNA,Letter RNA) where + mkBasepair (l,r) + | l >= A && l <= U && r >= A && r <= U + = BP $ 4 * getLetter l + getLetter r + | otherwise = NoBP + fromBasepair k + | k == NoBP || k == NS = (N,N) + | otherwise = let (l,r) = getBP k `divMod` 4 in (Letter l, Letter r) + {-# Inline mkBasepair #-} + {-# Inline fromBasepair #-} +-} + + +-- * Newtypes for extended secondary structures + +-- ** Encode which of three edges is engaged in base pairing + +-- | Each nucleotide in a pair may be paired using one of three edges: +-- watson-crick, sugar, or hoogsteen. + +newtype Edge = Edge {unEdge :: Int} + deriving (Eq,Ord,Ix,Generic) + +pattern W = Edge 0 +pattern S = Edge 1 +pattern H = Edge 2 + +instance Binary Edge +instance Serialize Edge +instance FromJSON Edge +instance ToJSON Edge + + + +-- | Human-readable Show instance. + +instance Show Edge where + show H = "H" + show S = "S" + show W = "W" + +-- | Human-readable Read instance. + +instance Read Edge where + readPrec = parens $ do + Ident s <- lexP + return $ case s of + "H" -> H + "S" -> S + "W" -> W + _ -> error $ "read Edge: " ++ s + +instance Bounded Edge where + minBound = W + maxBound = H + +instance Enum Edge where + toEnum = Edge + fromEnum = unEdge + +derivingUnbox "Edge" + [t| Edge -> Int |] [| unEdge |] [| Edge |] + +-- ** Is the base pair in cis or trans configuration + +-- | Nucleotides in a pairing may be in the cis(==?) or trans(==?) state. + +newtype CTisomerism = CT {unCT :: Int} + deriving (Eq,Ord,Ix,Generic) + +pattern Cis = CT 0 +pattern Trn = CT 1 + +instance Binary CTisomerism +instance Serialize CTisomerism +instance FromJSON CTisomerism +instance ToJSON CTisomerism + + +-- | Human-readable Show instance. + +instance Show CTisomerism where + show Cis = "C" + show Trn = "T" + +-- | Human-readable Read instance. + +instance Read CTisomerism where + readPrec = parens $ do + Ident s <- lexP + return $ case s of + "C" -> Cis + "T" -> Trn + _ -> error $ "read CTisomerism: " ++ s + +instance Bounded CTisomerism where + minBound = Cis + maxBound = Trn + +instance Enum CTisomerism where + toEnum = CT + fromEnum = unCT + +derivingUnbox "CTisomerism" + [t| CTisomerism -> Int |] [| unCT |] [| CT |] + + + +-- * Types + +-- | A basepair is simply a pair of Ints which are 0-indexing a sequence. + +type PairIdx = (Int,Int) + +-- | A pair as a tuple containing 'Nuc's. + +type Pair n = (Letter RNA n, Letter RNA n) + +-- | Annotation for a basepair. + +type ExtPairAnnotation = (CTisomerism,Edge,Edge) + +-- | An extended basepair is a basepair, annotated with edge and CTisomerism. + +type ExtPairIdx = (PairIdx,ExtPairAnnotation) + +-- | An extended basepair, with nucleotides an annotation. + +type ExtPair n = (Pair n, ExtPairAnnotation) + + + +-- * little helpers + +pattern CHH = (Cis,H,H) +pattern CHS = (Cis,H,S) +pattern CHW = (Cis,H,W) +pattern CSH = (Cis,S,H) +pattern CSS = (Cis,S,S) +pattern CSW = (Cis,S,W) +pattern CWH = (Cis,W,H) +pattern CWS = (Cis,W,S) +pattern CWW = (Cis,W,W) + +pattern THH = (Trn,H,H) +pattern THS = (Trn,H,S) +pattern THW = (Trn,H,W) +pattern TSH = (Trn,S,H) +pattern TSS = (Trn,S,S) +pattern TSW = (Trn,S,W) +pattern TWH = (Trn,W,H) +pattern TWS = (Trn,W,S) +pattern TWW = (Trn,W,W) + diff --git a/Biobase/Secondary/Constraint.hs b/Biobase/Secondary/Constraint.hs deleted file mode 100644 --- a/Biobase/Secondary/Constraint.hs +++ /dev/null @@ -1,115 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE StandaloneDeriving #-} - -module Biobase.Secondary.Constraint where - -import Data.Array.Repa.Index -import Data.Array.Repa.Shape -import Data.Char (toLower) -import Data.Primitive.Types -import qualified Data.Vector.Generic as VG -import qualified Data.Vector.Generic.Mutable as VGM -import qualified Data.Vector.Unboxed as VU - -import Data.PrimitiveArray -import Data.PrimitiveArray.Unboxed.Zero - -import Biobase.Secondary.Diagrams - - - --- | We can create a constraint from different sources - -class MkConstraint a where - mkConstraint :: a -> Constraint - --- | A constraint is nothing more than a vector of constraint characters --- together with a possible pairing for each character. - -newtype Constraint = Constraint {unConstraint :: VU.Vector (Char,Int)} - deriving (Show,Read,Eq) - -bonusCC = VU.fromList "()<>|" -{-# NOINLINE bonusCC #-} - -nobonusCC = VU.fromList ".x" -{-# NOINLINE nobonusCC #-} - --- | Given a 'Constraint', create an NxN matrix with bonus energies. These --- energies can be included in all pair-creating functions and will disallow or --- strongly favor certain pairings, while others will receive neither bonus nor --- malus. --- --- In case, a pair (i,j) is annotated as both, bonus- and malus-receiving, it --- will be set to receive a malus. This can happen, if something like "<" would --- give a bonus, but "x" gives a malus (and other cases). --- --- TODO and again, we should parametrize over "Energy", "Score", etc (that is, --- Prim a) - -bonusTable :: Double -> Double -> Constraint -> Arr0 DIM2 Double -bonusTable bonus malus (Constraint constraint) = arr where - arr = fromAssocs zeroDim (Z:.n:.n) 0 $ bonusBr ++ bonusAn ++ bonusBa ++ malusBr ++ malusAn ++ malusX - n = VU.length constraint -1 - infixl 1 `xor` - xor a b = a && not b || not a && b - -- "()" bonus energies - bonusBr = [ (Z:.i:.j,bonus) - | (i,('(',j)) <- zip [0..] $ VU.toList constraint - ] - malusBr = [ (Z:.i:.j,malus) - | i <- [0..n] - , j <- [i..n] - , let bi = constraint VU.! i - , let bj = constraint VU.! j - , fst bi == '(' && snd bi /= j || fst bj == ')' && snd bj /= i - ] - bonusAn = [ (Z:.i:.j,bonus) - | i<-[0..n] - , fst (constraint VU.! i) == '<' - , j<-[i+1..n] - ] ++ - [ (Z:.i:.j,bonus) - | j<-[0..n] - , fst (constraint VU.! j) == '>' - , i<-[0..j-1] - ] - malusAn = [ (Z:.i:.j,malus) - | i<-[0..n] - , j<-[i+1..n] - , fst (constraint VU.! j) == '<' - ] ++ - [ (Z:.i:.j,malus) - | i<-[0..n] - , j<-[i+1..n] - , fst (constraint VU.! i) == '>' - ] - bonusBa = [ (Z:.i:.j,bonus) - | i<-[0..n] - , j<-[i+1..n] - , fst (constraint VU.! i) == '|' || fst (constraint VU.! j) == '|' - ] - malusX = [ (Z:.i:.j,malus) - | i<-[0..n] - , j<-[i+1..n] - , fst (constraint VU.! i) == 'x' || fst (constraint VU.! j) == 'x' - ] - -{- -testC = putStrLn $ f as where - f [] = "" - f xs = show (take 9 xs) ++ "\n" ++ f (drop 9 xs) - as = toList $ bonusTable (1) 2 (mkConstraint "(<<..x|>)") --} - --- * Instances - -instance MkConstraint String where - mkConstraint xs = mkConstraint . VU.fromList . map toLower $ xs - -instance MkConstraint (VU.Vector Char) where - mkConstraint cs = Constraint $ VU.zip cs ks where - (D1S ks) = mkD1S cs diff --git a/Biobase/Secondary/Convert.hs b/Biobase/Secondary/Convert.hs new file mode 100644 --- /dev/null +++ b/Biobase/Secondary/Convert.hs @@ -0,0 +1,93 @@ + +-- | This module gives functionality to convert between different variants +-- of secondary structure elements. + +module Biobase.Secondary.Convert where + +import Biobase.Types.BioSequence + +import Biobase.Primary.Letter +import Biobase.Primary.Nuc.RNA +import Biobase.Secondary.Basepair +import Biobase.Secondary.Vienna (ViennaPair(..)) +import qualified Biobase.Secondary.Vienna as SV +import qualified Biobase.Secondary.Basepair as SB + + + +-- | @basepairConvert@ converts between different secondary structure base +-- pair representations. In general, the conversion is lossy, in particular +-- when "downsizing", say to @ViennaPair@. + +class BasepairConvert a b where + basepairConvert :: a -> b + + + +-- ** @(RNA,RNA) <-> Basepair@ + +instance BasepairConvert (Letter RNA n,Letter RNA n) Basepair where + basepairConvert (l,r) + | l >= A && l <= U && r >= A && r <= U + = BP $ 4 * getLetter l + getLetter r + | otherwise = NoBP + {-# Inline basepairConvert #-} + +instance BasepairConvert Basepair (Letter RNA n, Letter RNA n) where + basepairConvert k + | k == NoBP || k == NS = (N,N) + | otherwise = let (l,r) = getBP k `divMod` 4 in (Letter l, Letter r) + {-# Inline basepairConvert #-} + + + +-- ** @(RNA,RNA) <-> ViennaPair@ + +instance BasepairConvert (Letter RNA n, Letter RNA n) ViennaPair where + basepairConvert = \case + (C,G) -> SV.CG + (G,C) -> SV.GC + (G,U) -> SV.GU + (U,G) -> SV.UG + (A,U) -> SV.AU + (U,A) -> SV.UA + _ -> SV.NS + {-# Inline basepairConvert #-} + +instance BasepairConvert ViennaPair (Letter RNA n, Letter RNA n) where + basepairConvert = \case + SV.CG -> (C,G) + SV.GC -> (G,C) + SV.GU -> (G,U) + SV.UG -> (U,G) + SV.AU -> (A,U) + SV.UA -> (U,A) + SV.NS -> (N,N) + {-# Inline basepairConvert #-} + + + +-- ** @Basepair <-> ViennaPair@ + +instance BasepairConvert Basepair ViennaPair where + basepairConvert = \case + SB.AU -> SV.AU + SB.CG -> SV.CG + SB.GC -> SV.GC + SB.GU -> SV.GU + SB.UA -> SV.UA + SB.UG -> SV.UG + _ -> SV.NS + {-# Inline basepairConvert #-} + +instance BasepairConvert ViennaPair Basepair where + basepairConvert = \case + SV.AU -> SB.AU + SV.CG -> SB.CG + SV.GC -> SB.GC + SV.GU -> SB.GU + SV.UA -> SB.UA + SV.UG -> SB.UG + _ -> SB.NS + {-# Inline basepairConvert #-} + diff --git a/Biobase/Secondary/Diagrams.hs b/Biobase/Secondary/Diagrams.hs --- a/Biobase/Secondary/Diagrams.hs +++ b/Biobase/Secondary/Diagrams.hs @@ -1,51 +1,64 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE FlexibleInstances #-} -- | Types for RNA secondary structure. Types vary from the simplest array -- (D1Secondary) to rather complex ones. --- --- TODO The complex ones are still coming in from other libraries. --- --- TODO can we use Char8 instead of Char? --- --- TODO prepare for extended RNA secondary structures! -{-# LANGUAGE RecordWildCards #-} - module Biobase.Secondary.Diagrams where +import Control.Applicative +import Control.Arrow +import Control.Lens +import Data.Aeson +import Data.Binary +import Data.List ((\\)) +import Data.List (sort,groupBy,sortBy,intersperse) +import Data.List.Split (splitOn) +import Data.Serialize +import Data.Tuple.Select (sel1,sel2) +import Data.Tuple (swap) +import Data.Vector.Binary +import Data.Vector.Serialize +import GHC.Generics import qualified Data.Vector.Unboxed as VU -import Data.List (sort,groupBy,sortBy) -import Data.Tuple.Select (sel1,sel2) -import Data.Tuple (swap) -import Control.Arrow +import Text.Printf +import Control.DeepSeq -import Biobase.Primary -import Biobase.Secondary +import Biobase.Primary.Nuc +import Biobase.Secondary.Basepair -- | RNA secondary structure with 1-diagrams. Each nucleotide is paired with at --- most one other nucleotide. A nucleotide with index k \in [0..len-1] is --- paired if unD1S VU.! k > 0. +-- most one other nucleotide. A nucleotide with index @k@ in @[0..len-1]@ is +-- paired if @unD1S VU.! k >=0 0@ Unpaired status is @-1@. +-- +-- TODO Provide @iso@ between @D1Secondary@ and @RNAss@. newtype D1Secondary = D1S {unD1S :: VU.Vector Int} - deriving (Read,Show,Eq) + deriving (Read,Show,Eq,Generic,NFData) +instance Binary D1Secondary +instance Serialize D1Secondary +instance FromJSON D1Secondary +instance ToJSON D1Secondary + -- RNA secondary structure with 2-diagrams. Each nucleotide is paired with up -- to two other nucleotides. newtype D2Secondary = D2S {unD2S :: VU.Vector ( (Int,Edge,CTisomerism), (Int,Edge,CTisomerism) )} - deriving (Read,Show,Eq) + deriving (Read,Show,Eq,Generic) --- | +instance Binary D2Secondary +instance Serialize D2Secondary +instance FromJSON D2Secondary +instance ToJSON D2Secondary +-- | Conversion to and from 1-diagrams. + class MkD1Secondary a where mkD1S :: a -> D1Secondary fromD1S :: D1Secondary -> a --- | +-- | Conversion to and from 2-diagrams. class MkD2Secondary a where mkD2S :: a -> D2Secondary @@ -55,21 +68,21 @@ -- * Tree-based representation -- --- Tree -> d1/2Secondary ? +-- TODO Tree -> d1/2Secondary ? -- | A secondary-structure tree. Has no notion of pseudoknots. -data SSTree idx a = SSTree idx a [SSTree idx a] - | SSExt Int a [SSTree idx a] - deriving (Read,Show,Eq) +data SSTree idx a = SSTree idx a [SSTree idx a] + | SSExtern Int a [SSTree idx a] + deriving (Read,Show,Eq,Generic) -- | Create a tree from (pseudoknot-free [not checked]) 1-diagrams. d1sTree :: D1Secondary -> SSTree PairIdx () d1sTree s = ext $ sort ps where (len,ps) = fromD1S s - ext [] = SSExt len () [] - ext xs = SSExt len () . map tree $ groupBy (\l r -> snd l > fst r) xs -- ">=" would be partial allowance for 2-diagrams + ext [] = SSExtern len () [] + ext xs = SSExtern len () . map tree $ groupBy (\l r -> snd l > fst r) xs -- ">=" would be partial allowance for 2-diagrams tree [ij] = SSTree ij () [] tree (ij:xs) = SSTree ij () . map tree $ groupBy (\l r -> snd l > fst r) xs @@ -78,8 +91,8 @@ d2sTree :: D2Secondary -> SSTree ExtPairIdx () d2sTree s = ext $ sortBy d2Compare ps where (len,ps) = fromD2S s - ext [] = SSExt len () [] - ext xs = SSExt len () . map tree . groupBy d2Grouping $ xs + ext [] = SSExtern len () [] + ext xs = SSExtern len () . map tree . groupBy d2Grouping $ xs tree [ij] = SSTree ij () [] tree (ij:xs) = SSTree ij () . map tree . groupBy d2Grouping $ xs @@ -90,20 +103,6 @@ d2Grouping ((i,j),_) ((k,l),_) = i<=k && j>=l -test :: (Int,[ExtPairIdx]) -test = (20,test') - -test' = - [ ((2,15),(cis,wc,wc)) - , ((3,14),(cis,wc,wc)) - , ((4,13),(cis,wc,wc)) - , ((5,12),(cis,wc,wc)) - , ((6,10),(trans,wc,hoogsteen)) - , ((2,18),(trans,sugar,sugar)) - , ((15,18),(cis,sugar,sugar)) - ] - - -- * Instances for D1S -- | Conversion between D1S and D2S is lossy in D2S -> D1S @@ -126,7 +125,7 @@ -- TODO 'fromD2S' makes me wanna rewrite everything... instance MkD2Secondary D1Secondary where - mkD2S (D1S xs) = D2S . VU.map (\k -> ((k,wc,cis),(-1,unknownEdge,unknownCT))) $ xs + mkD2S = D2S . VU.map (\k -> ((k,W,Cis),(-1,W,Cis))) . unD1S fromD2S (D2S xs) = D1S . VU.map (sel1 . sel1) $ xs instance MkD2Secondary (Int,[ExtPairIdx]) where @@ -135,7 +134,7 @@ , (j, (i,e2,ct)) ]) ps f (x,y) z = if sel1 x == -1 then (z,y) else (x,z) - in D2S $ VU.accum f (VU.replicate len ((-1,unknownEdge,unknownCT),(-1,unknownEdge,unknownCT))) xs + in D2S $ VU.accum f (VU.replicate len ((-1,W,Cis),(-1,W,Cis))) xs fromD2S (D2S s) = ( VU.length s , let (xs,ys) = unzip . VU.toList $ s g i j = let z = s VU.! i in if sel1 (sel1 z) == j then sel2 (sel1 z) else sel2 (sel2 z) @@ -151,16 +150,12 @@ -- | A second primitive generator, requiring dictionary and String. This one -- generates pairs that are then used by the above instance. The dict is a list -- of possible brackets: ["()"] being the minimal set. --- --- NOTE no dictionary is returned by "fromD1S". --- --- TODO return dictionary that is actually seen? instance MkD1Secondary ([String],String) where mkD1S (dict,xs) = mkD1S (length xs,ps) where ps :: [(Int,Int)] - ps = dotBracket dict xs - fromD1S (D1S s) = ([], zipWith f [0..] $ VU.toList s) where + ps = unsafeDotBracket2pairlist dict xs + fromD1S (D1S s) = (["()"], zipWith f [0..] $ VU.toList s) where f k (-1) = '.' f k p | k>p = ')' @@ -175,22 +170,68 @@ -- | A "fast" instance for getting the pair list of vienna-structures. instance MkD1Secondary String where - mkD1S xs = mkD1S (["()"],xs) + mkD1S xs = mkD1S (["()" ::String],xs) fromD1S s = let (_::[String],res) = fromD1S s in res instance MkD1Secondary (VU.Vector Char) where - mkD1S xs = mkD1S (["()"],xs) + mkD1S xs = mkD1S (["()" ::String],xs) fromD1S s = let (_::[String],res::VU.Vector Char) = fromD1S s in res +-- * High-level parsing functionality for secondary structures + +-- | Completely canonical structure. +-- +-- TODO Check size of hairpins and interior loops? + +isCanonicalStructure :: String -> Bool +isCanonicalStructure = all (flip (elem @[]) "().") + +-- | Is constraint type structure, i.e. there can also be symbols present +-- that denote up- or downstream pairing. + +isConstraintStructure :: String -> Bool +isConstraintStructure = all (flip (elem @[]) "().<>{}|") + +-- | Take a structural string and split it into its constituents. +-- +-- If we decide to /NOT/ depend on @lens@ explicitly, another way to write +-- this is: +-- +-- @ +-- structures :: forall p f . (Profunctor p, Functor f) => p [String] (f [String]) -> p String (f String) +-- structures = dimap (splitOn "&") (fmap (concat . intersperse "&")) +-- @ + +structures :: Iso' String [String] +structures = iso (splitOn "&") (concat . intersperse "&") + +-- | A @fold@ structure is a single structure + +foldStructure :: Prism' String String +foldStructure = prism id to where + to s = case s^.structures of + [t] -> Right t + _ -> Left s + +-- | A @cofold@ structure has exactly two structures split by @&@ (which the +-- prism removes). + +cofoldStructure :: Prism' String (String,String) +cofoldStructure = prism from to where + from (l,r) = l ++ '&' : r + to s = case s^.structures of + [l,r] -> Right (l,r) + _ -> Left s + -- * Helper functions -- | Secondary structure parser which allows pseudoknots, if they use different -- kinds of brackets. -dotBracket :: [String] -> String -> [(Int,Int)] -dotBracket dict xs = sort . concatMap (f xs) $ dict where +unsafeDotBracket2pairlist :: [String] -> String -> [(Int,Int)] +unsafeDotBracket2pairlist dict xs = sort . concatMap (f xs) $ dict where f xs [l,r] = g 0 [] . map (\x -> if x `elem` [l,r] then x else '.') $ xs where g :: Int -> [Int] -> String -> [(Int,Int)] g _ st [] = [] @@ -200,3 +241,54 @@ g k (s:st) (x:xs) | r==x = (s,k) : g (k+1) st xs g a b c = error $ show (a,b,c) + +-- | Secondary structure parser with a notion of errors. We either return a +-- @Right@ structure, including flags, or a @Left@ error. + +dotBracket2pairlist :: [String] -> String -> Either String ( [(Int,Int)] ) +dotBracket2pairlist dict str = fmap (sort . concat) . sequence . map (f str) $ dict where + f ys [l,r] = g 0 [] . map (\x -> if x `elem` [l,r] then x else '.') $ ys where + g :: Int -> [Int] -> String -> Either String ( [(Int,Int)] ) + g _ [] [] = pure [] + g k st ('.':xs) = g (k+1) st xs + g k st (x:xs) | l==x = g (k+1) (k:st) xs + g k (s:st) (x:xs) | r==x = ((s,k):) <$> g (k+1) st xs + g k [] xs = Left $ printf "too many closing brackets at position %d: '%s' (dot-bracket: %s)" k xs str + g k st [] = Left $ printf "too many opening brackets, opening bracket(s) at: %s (dot-bracket: %s)" (show $ reverse st) str + g a b c = Left $ printf "unspecified error: %s (dot-bracket: %s)" (show (a,b,c)) str + f xs lr@(_:_:_:_) = Left $ printf "unsound dictionary: %s (dot-bracket: %s)" lr str + f xs lr = Left $ printf "unspecified error: dict: %s, input: %s (dot-bracket: %s)" lr xs str + +-- | Calculates the distance between two vienna strings. + +viennaStringDistance :: Bool -> Bool -> String -> String -> (String,Int) +viennaStringDistance sPairs tPairs s t = (t,length $ ss++tt) where + s' = either error id . dotBracket2pairlist ["()"] $ s + t' = either error id . dotBracket2pairlist ["()"] $ t + ss = if sPairs then s' \\ t' else [] + tt = if tPairs then t' \\ s' else [] + +-- | Calculate the distance between two 'D1Secondary' structures, that live +-- in the same underlying space. In particular, this probably only works +-- for structures on the same primary sequence. +-- +-- This function assumes somewhat dense structures, as it is @O(2n)@ with +-- @n@ the length of the underlying vectors. +-- +-- @(i,k)@ vs @(j,l)@ +-- +-- TODO error out on weird inputs! + +d1Distance :: D1Secondary -> D1Secondary -> Int +d1Distance (D1S x) (D1S y) +-- | VU.length x /= VU.length y = error "d1Distance called on vectors with differing lengths!" + | otherwise = (`div` 2) . VU.sum $ VU.zipWith chk (x VU.++ xx) (y VU.++ yy) + where xx = VU.replicate (VU.length y - VU.length x) (-2) + yy = VU.replicate (VU.length x - VU.length y) (-2) + chk i j | i==j = 0 + | i < 0 && j < 0 = 0 + | i >= 0 && j >= 0 = 2 + | otherwise = 1 + {-# Inline chk #-} +{-# NoInline d1Distance #-} + diff --git a/Biobase/Secondary/Isostericity.hs b/Biobase/Secondary/Isostericity.hs --- a/Biobase/Secondary/Isostericity.hs +++ b/Biobase/Secondary/Isostericity.hs @@ -1,28 +1,28 @@ -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE TemplateHaskell #-} -- | Provides detailed information on isostericity of RNA basepairs. All data -- is extracted from csv files which were created from supplemental files in: -- +-- @ -- Frequency and isostericity of RNA base pairs -- Jesse Stombaugh, Craig L. Zirbel, Eric Westhof, and Neocles B. Leontis -- Nucl. Acids Res. (2009) -- doi:10.1093/nar/gkp011 +-- @ +-- module Biobase.Secondary.Isostericity where -import Data.ByteString.Char8 (ByteString) -import Data.FileEmbed (embedFile) -import Data.Function (on) -import Data.List +import Data.ByteString.Char8 (ByteString) +import Data.FileEmbed (makeRelativeToProject, embedFile) +import Data.Function (on) +import Data.List +import Data.Tuple.Select import qualified Data.ByteString.Char8 as BS import qualified Data.Map as M -import Text.CSV +import Text.CSV -import Biobase.Primary -import Biobase.Secondary +import Biobase.Primary.Nuc +import Biobase.Secondary.Basepair @@ -43,7 +43,7 @@ -- -- TODO inClass missing -instance IsostericityLookup ExtPair where +instance IsostericityLookup (ExtPair n) where getClasses p | Just cs <- M.lookup p defaultIsostericityMap = cs @@ -54,18 +54,17 @@ -- -- TODO inClass missing -instance IsostericityLookup Pair where +instance IsostericityLookup (Pair n) where getClasses p - | Just cs <- M.lookup (p,cWW) defaultIsostericityMap + | Just cs <- M.lookup (p,CWW) defaultIsostericityMap = cs | otherwise = [] - inClass x = map (baseP.fst) -- remove extended information - . filter ((cWW==).baseT.fst) -- keep only cWW pairs (baseT-ype) - . filter ((x `elem`).snd) -- select based on class + inClass x = map (sel1 . fst) -- remove extended information + . filter ((CWW==). snd . fst) -- keep only cWW pairs (baseT-ype) + . filter ((x `elem`).snd) -- select based on class $ M.assocs defaultIsostericityMap - -- ** default data -- | The default isostericity mapping. @@ -78,20 +77,23 @@ -- | Process CSV list-of-lists to get the isostericity data. -mkIsostericityList :: [[[String]]] -> [(ExtPair, [String])] +mkIsostericityList :: [[[String]]] -> [(ExtPair n, [String])] mkIsostericityList gs = nubBy ((==) `on` fst) . concatMap turn . concatMap f $ gs where f g = map (\e -> ( ( let [x,y] = fst e - in (mkNuc x, mkNuc y), threeChar bpt + in (charRNA x, charRNA y), read bpt ) , nub $ snd e) ) $ map entry xs where bpt = head $ head g xs = tail g - entry x = (x!!0, map (filter (\z -> not $ z `elem` "()")) . takeWhile ('I' `elem`) . drop 2 $ x) + entry x = (x!!0, map (filter (\z -> not $ z `elem` bracket)) . takeWhile ('I' `elem`) . drop 2 $ x) + bracket :: String + bracket = "()" turn entry@(((x,y),(wc,tx,ty)), cs) = [entry, (((y,x),(wc,ty,tx)), cs)] -- | Simple parsing of raw CSV data. +parsedCSV :: [[[Field]]] parsedCSV = filter (not . null) gs where gs = map (filter ((""/=).head)) . groupBy (\x y -> ""/= (head y)) $ csv Right csv = parseCSV "isostericity/detailed" $ BS.unpack detailedCSV @@ -103,5 +105,5 @@ -- | Raw CSV data, embedded into the library. detailedCSV :: ByteString -detailedCSV = $(embedFile "sources/isostericity-detailed.csv") +detailedCSV = $(makeRelativeToProject "sources/isostericity-detailed.csv" >>= embedFile) diff --git a/Biobase/Secondary/New.hs b/Biobase/Secondary/New.hs new file mode 100644 --- /dev/null +++ b/Biobase/Secondary/New.hs @@ -0,0 +1,101 @@ + +-- | New parsers and structures for secondary structures. The structures here a strict. +-- +-- TODO Parser should check if a @#Vienna Secondary Structure@ or @#Extended Secondary Structure@ precedes the entries. + +module Biobase.Secondary.New where + +import Control.Applicative +import Control.Lens +import Control.Monad.Except +import Data.Attoparsec.ByteString.Char8 +import Data.ByteString.Char8 (ByteString,pack) +import Data.Functor +import Data.Tree +import Data.Vector (Vector, fromList) +import GHC.Generics (Generic) + + + +-- | A completely closed sub-structure. An unpaired region @.@ is closed. A +-- paired region @(r)@ is closed, where @r@ contains arbitrarily many unpaired +-- and paired elements. +-- +-- TODO Should be extended with @Extended@, but this requires knowing which of +-- the ends overlap with paired: left, right, or both. + +data SubStructure (t :: *) a + = Unpaired { _label :: !a } + | Paired { _label :: !a, _subStructures :: !(Vector (SubStructure t a)) } + deriving (Show, Read, Functor, Traversable, Foldable, Generic, Eq, Ord) +makeLenses ''SubStructure +makePrisms ''SubStructure + +-- | A full structure is composed of a number of sub-structures. The empty +-- structure is a full structure. + +newtype FullStructure (t :: *) a + = FullStructure { _fullStructure :: Vector (SubStructure t a) } + deriving (Show, Read, Functor, Traversable, Foldable, Generic, Eq, Ord) +makeLenses ''FullStructure + + + +-- ** Parses a ViennaRNA secondary structure string. + +pUnpaired :: Parser (SubStructure () ()) +pUnpaired = Unpaired () <$ char '.' +{-# Inlinable pUnpaired #-} + +pPaired :: Parser (SubStructure () ()) +pPaired = Paired () <$ char '(' <*> (fromList <$> many pSubStructure) <* char ')' +{-# Inlinable pPaired #-} + +pSubStructure :: Parser (SubStructure () ()) +pSubStructure = pUnpaired <|> pPaired +{-# Inlinable pSubStructure #-} + +pFullStructure :: Parser (FullStructure () ()) +pFullStructure = FullStructure <$> fromList <$> many pSubStructure <* endOfInput +{-# Inlinable pFullStructure #-} + +newtype StructureParseError = StructureParseError String + deriving (Show) + +parseVienna :: MonadError StructureParseError m ⇒ ByteString -> m (FullStructure () ()) +parseVienna = either (throwError . StructureParseError) return . parseOnly pFullStructure +{-# Inlinable parseVienna #-} + + + +-- ** Transform into a @Tree@. + +-- | Transform a 'FullStructure' into a 'Tree'. +-- +-- Given a full structure generated like this: +-- @ +-- s = either (error . show) id $ parseVienna $ pack ".()(())." +-- @ +-- +-- a tree of just the base paired can be created with +-- @ +-- toTree (preview (_Paired._1)) () s +-- @ + +toTree + :: (SubStructure t a -> Maybe b) + -- ^ how to handle substructure elements? @Nothing@ means discard this + -- substructure and all children. + -> b + -- ^ The root label + -> FullStructure (t :: *) a + -- ^ The @FullStructure@ to transform into a @Tree@. + -> Tree b +toTree f r (FullStructure ts) = Node r $ fmap go ts ^.. traverse . _Just + where + go u@Unpaired{} = (`Node` []) <$> f u + go p@Paired{} = case f p of + Nothing -> Nothing + Just lbl -> Just $ Node lbl $ (fmap go $ p^.subStructures) ^.. traverse . _Just +{-# Inlinable toTree #-} + diff --git a/Biobase/Secondary/PseudoKnots.hs b/Biobase/Secondary/PseudoKnots.hs deleted file mode 100644 --- a/Biobase/Secondary/PseudoKnots.hs +++ /dev/null @@ -1,61 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE BangPatterns #-} - --- | Methods to transform a secondary structure containing pseudoknots into a --- structure which is pseudoknot-free. --- --- TODO Until a better name is found, this module is home to functions for --- "de-pseudoknotting" structures. --- --- TODO Check if there are corner-cases remaining when considering 2-diagrams. - -module Biobase.Secondary.PseudoKnots where - -import qualified Data.Vector.Unboxed as VU -import Data.List - -import Biobase.Secondary - - - --- | Try to removed pseudoknots from the "pairlist". This works by counting for --- each pair, how many pairs are incompatible with it. Then those with most --- incompatibilities are successively removed. This function might well remove --- more than necessary! - -class RemovePseudoKnots a where - removeByCounting :: a -> a - --- | Remove pseudoknotted pairs from RNA secondary structures. - -instance RemovePseudoKnots (VU.Vector PairIdx) where - removeByCounting = VU.force . wrapRemove where - wrapRemove !ps - | VU.null cnts = ps -- there are no pairs - | mmx == 0 = ps -- there are no incompatibilities - | otherwise = wrapRemove $ VU.take pos ps VU.++ VU.drop (pos+1) ps - where - cnts = VU.map incomp ps - mmx = VU.maximum cnts - Just pos = VU.elemIndex mmx cnts - incomp (i,j) = VU.length $ VU.filter (\(k,l) -> i i a + +-- | Remove pseudoknotted pairs from RNA secondary structures. + +instance RemovePseudoKnots (VU.Vector PairIdx) where + removeByCounting = VU.force . wrapRemove where + wrapRemove !ps + | VU.null cnts = ps -- there are no pairs + | mmx == 0 = ps -- there are no incompatibilities + | otherwise = wrapRemove $ VU.take pos ps VU.++ VU.drop (pos+1) ps + where + cnts = VU.map incomp ps + mmx = VU.maximum cnts + Just pos = VU.elemIndex mmx cnts + incomp (i,j) = VU.length $ VU.filter (\(k,l) -> i i Shape (sh :. ViennaPair) where - rank (sh:._) = rank sh + 1 - zeroDim = zeroDim:.ViennaPair 0 - unitDim = unitDim:.ViennaPair 1 -- TODO does this one make sense? - intersectDim (sh1:.n1) (sh2:.n2) = intersectDim sh1 sh2 :. min n1 n2 - addDim (sh1:.ViennaPair n1) (sh2:.ViennaPair n2) = addDim sh1 sh2 :. ViennaPair (n1+n2) -- TODO will not necessarily yield a valid ViennaPair - size (sh1:.ViennaPair n) = size sh1 * n - sizeIsValid (sh1:.ViennaPair n) = sizeIsValid (sh1:.n) - toIndex (sh1:.ViennaPair sh2) (sh1':.ViennaPair sh2') = toIndex (sh1:.sh2) (sh1':.sh2') - fromIndex (ds:.ViennaPair d) n = fromIndex ds (n `quotInt` d) :. ViennaPair r where - r | rank ds == 0 = n - | otherwise = n `remInt` d - inShapeRange (sh1:.n1) (sh2:.n2) (idx:.i) = i>=n1 && i error "empty list in shapeOfList/Primary" - x:xs -> shapeOfList xs :. ViennaPair x - deepSeq (sh:.n) x = deepSeq sh (n `seq` x) - {-# INLINE rank #-} - {-# INLINE zeroDim #-} - {-# INLINE unitDim #-} - {-# INLINE intersectDim #-} - {-# INLINE addDim #-} - {-# INLINE size #-} - {-# INLINE sizeIsValid #-} - {-# INLINE toIndex #-} - {-# INLINE fromIndex #-} - {-# INLINE inShapeRange #-} - {-# INLINE listOfShape #-} - {-# INLINE shapeOfList #-} - {-# INLINE deepSeq #-} +derivingUnbox "ViennaPair" + [t| ViennaPair -> Int |] + [| unViennaPair |] + [| ViennaPair |] -instance (Eq sh, Shape sh, Show sh, ExtShape sh) => ExtShape (sh :. ViennaPair) where - subDim (sh1:.ViennaPair n1) (sh2:.ViennaPair n2) = subDim sh1 sh2 :. (ViennaPair $ n1-n2) - rangeList (sh1:.ViennaPair n1) (sh2:.ViennaPair n2) = [sh:.ViennaPair n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2)]] +instance Binary (ViennaPair) +instance Serialize (ViennaPair) +instance FromJSON (ViennaPair) +instance ToJSON (ViennaPair) -(vpNP:vpCG:vpGC:vpGU:vpUG:vpAU:vpUA:vpNS:vpUndefined:_) = map ViennaPair [0..] +instance Index ViennaPair where + data LimitType ViennaPair + = Canonical | Extended + linearIndex _ (ViennaPair p) = p + {-# Inline linearIndex #-} + size h = case h of { Canonical → 7; Extended → 9 } + {-# Inline size #-} + inBounds h (ViennaPair p) = 0 <= p && p < size h + {-# Inline inBounds #-} +instance IndexStream z => IndexStream (z:.ViennaPair) where + streamUp (ls:..l) (hs:..h) = flatten mk step $ streamUp ls hs + where mk z = return (z,size l - 1) + step (z,k) + | k > size h -1 = return $ Done + | otherwise = return $ Yield (z:.ViennaPair k) (z,k+1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamUp #-} + streamDown (ls:..l) (hs:..h) = flatten mk step $ streamDown ls hs + where mk z = return (z,size h - 1) + step (z,k) + | k < size l -1 = return $ Done + | otherwise = return $ Yield (z:.ViennaPair k) (z,k-1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamDown #-} + +instance IndexStream ViennaPair where + + + +pattern NP = ViennaPair 0 :: ViennaPair +pattern CG = ViennaPair 1 :: ViennaPair +pattern GC = ViennaPair 2 :: ViennaPair +pattern GU = ViennaPair 3 :: ViennaPair +pattern UG = ViennaPair 4 :: ViennaPair +pattern AU = ViennaPair 5 :: ViennaPair +pattern UA = ViennaPair 6 :: ViennaPair +-- | Non-standard base pair +pattern NS = ViennaPair 7 :: ViennaPair +pattern Undef = ViennaPair 8 :: ViennaPair + +{- class MkViennaPair a where mkViennaPair :: a -> ViennaPair fromViennaPair :: ViennaPair -> a -instance MkViennaPair (Nuc,Nuc) where - mkViennaPair (b1,b2) = viennaPairTable `PA.index` (Z:.b1:.b2) - {- - | b1==nC&&b2==nG = vpCG - | b1==nG&&b2==nC = vpGC - | b1==nG&&b2==nU = vpGU - | b1==nU&&b2==nG = vpUG - | b1==nA&&b2==nU = vpAU - | b1==nU&&b2==nA = vpUA - | otherwise = vpNS - -} +instance MkViennaPair (Letter RNA, Letter RNA) where + mkViennaPair = \case + (C,G) -> CG + (G,C) -> GC + (G,U) -> GU + (U,G) -> UG + (A,U) -> AU + (U,A) -> UA + _ -> NS {-# INLINE mkViennaPair #-} - fromViennaPair p - | p==vpCG = (nC,nG) - | p==vpGC = (nG,nC) - | p==vpGU = (nG,nU) - | p==vpUG = (nU,nG) - | p==vpAU = (nA,nU) - | p==vpUA = (nU,nA) - | otherwise = error "non-standard pairs can't be backcasted" + fromViennaPair = \case + CG -> (C,G) + GC -> (G,C) + GU -> (G,U) + UG -> (U,G) + AU -> (A,U) + UA -> (U,A) + _ -> error "non-standard pairs can't be backcasted" {-# INLINE fromViennaPair #-} +-} -viennaPairTable :: Arr0 (Z:.Nuc:.Nuc) ViennaPair -viennaPairTable = fromAssocs (Z:.nN:.nN) (Z:.nU:.nU) vpNS - [ (Z:.nC:.nG , vpCG) - , (Z:.nG:.nC , vpGC) - , (Z:.nG:.nU , vpGU) - , (Z:.nU:.nG , vpUG) - , (Z:.nA:.nU , vpAU) - , (Z:.nU:.nA , vpUA) +isViennaPair :: Letter RNA m -> Letter RNA n -> Bool +isViennaPair l r = l==C && r==G + || l==G && r==C + || l==A && r==U + || l==U && r==A + || l==G && r==U + || l==U && r==G +{-# INLINE isViennaPair #-} + +viennaPairTable :: Unboxed (Z:.Letter RNA n:.Letter RNA n) ViennaPair +viennaPairTable = fromAssocs (ZZ:..LtLetter maxBound:..LtLetter maxBound) NS + [ (Z:.C:.G , CG) + , (Z:.G:.C , GC) + , (Z:.G:.U , GU) + , (Z:.U:.G , UG) + , (Z:.A:.U , AU) + , (Z:.U:.A , UA) ] {-# NOINLINE viennaPairTable #-} -deriving instance VGM.MVector VU.MVector ViennaPair -deriving instance VG.Vector VU.Vector ViennaPair -deriving instance VU.Unbox ViennaPair -deriving instance Prim ViennaPair - instance Enum ViennaPair where toEnum x | x>=0 && x<=7 = ViennaPair x @@ -121,14 +141,8 @@ {-# INLINE fromEnum #-} instance Bounded ViennaPair where - minBound = vpNP - maxBound = vpNS - -instance Bounds ViennaPair where - minNormal = vpCG - maxNormal = vpUA - minExtended = vpNP - maxExtended = vpNS + minBound = NP + maxBound = NS instance Show ViennaPair where show x @@ -142,27 +156,28 @@ | x ==' ' = readsPrec p (y:xs) | Just n <- (x:y:[]) `lookup` s2p = [(n,xs)] | otherwise = [] - where s2p = (map swap pairToString) + where s2p = (P.map swap pairToString) -- | reverse a vienna pair revPair :: ViennaPair -> ViennaPair -revPair p - | p==vpCG = vpGC - | p==vpGC = vpCG - | p==vpGU = vpUG - | p==vpUG = vpGU - | p==vpAU = vpUA - | p==vpUA = vpAU - | p==vpNP = vpNP - | p==vpNS = vpNS +revPair = \case + CG -> GC + GC -> CG + GU -> UG + UG -> GU + AU -> UA + UA -> AU + NP -> NP + NS -> NS -- * Convenience structures -cguaP = [vpCG..vpUA] -cgnsP = [vpCG..vpNS] -pairToString = [(vpCG,"CG"),(vpGC,"GC"),(vpUA,"UA"),(vpAU,"AU"),(vpGU,"GU"),(vpUG,"UG"),(vpNS,"NS"),(vpNP,"NP")] +cguaP = [CG .. UA] +cgnsP = [CG .. NS] +pairToString = [(CG,"CG"),(GC,"GC"),(UA,"UA"),(AU,"AU"),(GU,"GU"),(UG,"UG"),(NS,"NS"),(NP,"NP")] + diff --git a/BiobaseXNA.cabal b/BiobaseXNA.cabal --- a/BiobaseXNA.cabal +++ b/BiobaseXNA.cabal @@ -1,69 +1,184 @@ +cabal-version: 2.2 name: BiobaseXNA -version: 0.6.2.0 +version: 0.11.1.1 author: Christian Hoener zu Siederdissen -maintainer: choener@tbi.univie.ac.at -homepage: http://www.tbi.univie.ac.at/~choener/ -copyright: Christian Hoener zu Siederdissen, 2011-2012 +maintainer: choener@bioinf.uni-leipzig.de +homepage: https://github.com/choener/BiobaseXNA +bug-reports: https://github.com/choener/BiobaseXNA/issues +copyright: Christian Hoener zu Siederdissen, 2011 - 2021 category: Bioinformatics -synopsis: Efficient RNA/DNA representations -license: GPL-3 +synopsis: Efficient RNA/DNA/Protein Primary/Secondary Structure +license: BSD-3-Clause license-file: LICENSE build-type: Simple stability: experimental -cabal-version: >= 1.6.0 +tested-with: GHC == 8.8, GHC == 8.10, GHC == 9.0 description: This is a base library for bioinformatics with emphasis on RNA - and DNA primary structure and related tools. Provided are - efficient encodings for short sequences, as required by RNA - folding tools. Extended RNA secondary structures can be - represented as well. + and DNA primary structure as well as amino acid sequences. . + Provided are efficient encodings for short (limited by the + amount of RAM) sequences, as required by RNA folding tools. + Extended RNA secondary structures can be represented as well. + . Contains data from: . + @ Frequency and isostericity of RNA base pairs - . Jesse Stombaugh, Craig L. Zirbel, Eric Westhof, and Neocles B. Leontis - . Nucl. Acids Res. (2009) + @ . - . - . - . - New in 0.6.2.0 - . - * Updated to PrimitiveArray >= 0.2.0.0 + + extra-source-files: sources/isostericity-matrices.csv sources/isostericity-detailed.csv + sources/iupac-nucleotides + sources/codontable + changelog.md + README.md -library - build-depends: - base >3 && <5, - containers, - bytestring, - csv, - file-embed, - primitive, - text, - tuple, - vector >=0.9 && <0.10, - PrimitiveArray == 0.2.0.0 +data-files: + sources/iupac-nucleotides + +common deps + build-depends: base >= 4.7 && < 5.0 + , aeson >= 1.0 + , attoparsec >= 0.13 + , binary >= 0.7 + , bytes >= 0.15 + , bytestring >= 0.10 + , cereal >= 0.4 + , cereal-vector >= 0.2 + , containers >= 0.5 + , csv >= 0.1 + , data-default >= 0.7 + , deepseq >= 1.3 + , file-embed >= 0.0.8 + , hashable >= 1.2 + , lens >= 4.0 + , mtl >= 2.0 + , primitive >= 0.5 + , QuickCheck >= 2.7 + , split >= 0.2 + , text >= 1.0 + , tuple >= 0.3 + , vector >= 0.11 + , vector-binary-instances >= 0.2 + , vector-th-unbox >= 0.2 + -- + , bimaps == 0.1.0.* + , BiobaseENA == 0.0.0.* + , BiobaseTypes == 0.2.1.* + , DPutils == 0.1.1.* + , ForestStructures == 0.0.1.* + , PrimitiveArray == 0.10.1.* + default-extensions: BangPatterns + , DataKinds + , DeriveDataTypeable + , DeriveFunctor + , DeriveGeneric + , DeriveGeneric + , DeriveTraversable + , EmptyDataDecls + , FlexibleContexts + , FlexibleInstances + , GeneralizedNewtypeDeriving + , LambdaCase + , MultiParamTypeClasses + , PatternSynonyms + , PolyKinds + , RankNTypes + , RecordWildCards + , ScopedTypeVariables + , StandaloneDeriving + , TemplateHaskell + , TypeApplications + , TypeFamilies + , TypeOperators + , UndecidableInstances + , UnicodeSyntax + , ViewPatterns + default-language: + Haskell2010 + ghc-options: + -O2 -funbox-strict-fields + + +library + import: + deps exposed-modules: + -- new + Biobase.Primary.AA + Biobase.Primary.Letter + Biobase.Primary.Nuc.RNA + Biobase.Primary.Unknown + -- old Biobase.Primary Biobase.Primary.Bounds Biobase.Primary.Hashed + Biobase.Primary.IUPAC + Biobase.Primary.Nuc + Biobase.Primary.Nuc.Conversion + Biobase.Primary.Nuc.DNA + Biobase.Primary.Nuc.XNA + Biobase.Primary.Pretty + Biobase.Primary.Trans Biobase.Secondary - Biobase.Secondary.Constraint + Biobase.Secondary.New + Biobase.Secondary.Basepair +-- Biobase.Secondary.Constraint + Biobase.Secondary.Convert Biobase.Secondary.Diagrams Biobase.Secondary.Isostericity - Biobase.Secondary.PseudoKnots + Biobase.Secondary.Pseudoknots + Biobase.Secondary.Structure Biobase.Secondary.Vienna + + +executable SubOptDistance + import: + deps + build-depends: base + , BiobaseXNA + , cmdargs >= 0.10 + main-is: + SubOptDistance.hs + hs-source-dirs: + src + + + +test-suite properties + import: + deps + type: + exitcode-stdio-1.0 + main-is: + properties.hs ghc-options: - -Odph -funbox-strict-fields -fspec-constr + -threaded -rtsopts -with-rtsopts=-N + hs-source-dirs: + tests + default-language: + Haskell2010 + default-extensions: TemplateHaskell + , ScopedTypeVariables + build-depends: base + , QuickCheck >= 2.7 + , tasty >= 0.11 + , tasty-quickcheck >= 0.8 + , tasty-th >= 0.1 + -- + , BiobaseXNA + + source-repository head type: git diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,675 +1,30 @@ - GNU GENERAL PUBLIC LICENSE - Version 3, 29 June 2007 - - Copyright (C) 2007 Free Software Foundation, Inc. - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The GNU General Public License is a free, copyleft license for -software and other kinds of works. - - The licenses for most software and other practical works are designed -to take away your freedom to share and change the works. By contrast, -the GNU General Public License is intended to guarantee your freedom to -share and change all versions of a program--to make sure it remains free -software for all its users. We, the Free Software Foundation, use the -GNU General Public License for most of our software; it applies also to -any other work released this way by its authors. You can apply it to -your programs, too. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -them if you wish), that you receive source code or can get it if you -want it, that you can change the software or use pieces of it in new -free programs, and that you know you can do these things. - - To protect your rights, we need to prevent others from denying you -these rights or asking you to surrender the rights. Therefore, you have -certain responsibilities if you distribute copies of the software, or if -you modify it: responsibilities to respect the freedom of others. - - For example, if you distribute copies of such a program, whether -gratis or for a fee, you must pass on to the recipients the same -freedoms that you received. You must make sure that they, too, receive -or can get the source code. And you must show them these terms so they -know their rights. - - Developers that use the GNU GPL protect your rights with two steps: -(1) assert copyright on the software, and (2) offer you this License -giving you legal permission to copy, distribute and/or modify it. - - For the developers' and authors' protection, the GPL clearly explains -that there is no warranty for this free software. For both users' and -authors' sake, the GPL requires that modified versions be marked as -changed, so that their problems will not be attributed erroneously to -authors of previous versions. - - Some devices are designed to deny users access to install or run -modified versions of the software inside them, although the manufacturer -can do so. This is fundamentally incompatible with the aim of -protecting users' freedom to change the software. The systematic -pattern of such abuse occurs in the area of products for individuals to -use, which is precisely where it is most unacceptable. Therefore, we -have designed this version of the GPL to prohibit the practice for those -products. If such problems arise substantially in other domains, we -stand ready to extend this provision to those domains in future versions -of the GPL, as needed to protect the freedom of users. - - Finally, every program is threatened constantly by software patents. -States should not allow patents to restrict development and use of -software on general-purpose computers, but in those that do, we wish to -avoid the special danger that patents applied to a free program could -make it effectively proprietary. To prevent this, the GPL assures that -patents cannot be used to render the program non-free. - - The precise terms and conditions for copying, distribution and -modification follow. - - TERMS AND CONDITIONS - - 0. Definitions. - - "This License" refers to version 3 of the GNU General Public License. - - "Copyright" also means copyright-like laws that apply to other kinds of -works, such as semiconductor masks. - - "The Program" refers to any copyrightable work licensed under this -License. Each licensee is addressed as "you". "Licensees" and -"recipients" may be individuals or organizations. - - To "modify" a work means to copy from or adapt all or part of the work -in a fashion requiring copyright permission, other than the making of an -exact copy. The resulting work is called a "modified version" of the -earlier work or a work "based on" the earlier work. - - A "covered work" means either the unmodified Program or a work based -on the Program. - - To "propagate" a work means to do anything with it that, without -permission, would make you directly or secondarily liable for -infringement under applicable copyright law, except executing it on a -computer or modifying a private copy. Propagation includes copying, -distribution (with or without modification), making available to the -public, and in some countries other activities as well. - - To "convey" a work means any kind of propagation that enables other -parties to make or receive copies. Mere interaction with a user through -a computer network, with no transfer of a copy, is not conveying. - - An interactive user interface displays "Appropriate Legal Notices" -to the extent that it includes a convenient and prominently visible -feature that (1) displays an appropriate copyright notice, and (2) -tells the user that there is no warranty for the work (except to the -extent that warranties are provided), that licensees may convey the -work under this License, and how to view a copy of this License. If -the interface presents a list of user commands or options, such as a -menu, a prominent item in the list meets this criterion. - - 1. Source Code. - - The "source code" for a work means the preferred form of the work -for making modifications to it. "Object code" means any non-source -form of a work. - - A "Standard Interface" means an interface that either is an official -standard defined by a recognized standards body, or, in the case of -interfaces specified for a particular programming language, one that -is widely used among developers working in that language. - - The "System Libraries" of an executable work include anything, other -than the work as a whole, that (a) is included in the normal form of -packaging a Major Component, but which is not part of that Major -Component, and (b) serves only to enable use of the work with that -Major Component, or to implement a Standard Interface for which an -implementation is available to the public in source code form. A -"Major Component", in this context, means a major essential component -(kernel, window system, and so on) of the specific operating system -(if any) on which the executable work runs, or a compiler used to -produce the work, or an object code interpreter used to run it. - - The "Corresponding Source" for a work in object code form means all -the source code needed to generate, install, and (for an executable -work) run the object code and to modify the work, including scripts to -control those activities. However, it does not include the work's -System Libraries, or general-purpose tools or generally available free -programs which are used unmodified in performing those activities but -which are not part of the work. For example, Corresponding Source -includes interface definition files associated with source files for -the work, and the source code for shared libraries and dynamically -linked subprograms that the work is specifically designed to require, -such as by intimate data communication or control flow between those -subprograms and other parts of the work. - - The Corresponding Source need not include anything that users -can regenerate automatically from other parts of the Corresponding -Source. - - The Corresponding Source for a work in source code form is that -same work. - - 2. Basic Permissions. - - All rights granted under this License are granted for the term of -copyright on the Program, and are irrevocable provided the stated -conditions are met. This License explicitly affirms your unlimited -permission to run the unmodified Program. The output from running a -covered work is covered by this License only if the output, given its -content, constitutes a covered work. This License acknowledges your -rights of fair use or other equivalent, as provided by copyright law. - - You may make, run and propagate covered works that you do not -convey, without conditions so long as your license otherwise remains -in force. You may convey covered works to others for the sole purpose -of having them make modifications exclusively for you, or provide you -with facilities for running those works, provided that you comply with -the terms of this License in conveying all material for which you do -not control copyright. Those thus making or running the covered works -for you must do so exclusively on your behalf, under your direction -and control, on terms that prohibit them from making any copies of -your copyrighted material outside their relationship with you. - - Conveying under any other circumstances is permitted solely under -the conditions stated below. Sublicensing is not allowed; section 10 -makes it unnecessary. - - 3. Protecting Users' Legal Rights From Anti-Circumvention Law. - - No covered work shall be deemed part of an effective technological -measure under any applicable law fulfilling obligations under article -11 of the WIPO copyright treaty adopted on 20 December 1996, or -similar laws prohibiting or restricting circumvention of such -measures. - - When you convey a covered work, you waive any legal power to forbid -circumvention of technological measures to the extent such circumvention -is effected by exercising rights under this License with respect to -the covered work, and you disclaim any intention to limit operation or -modification of the work as a means of enforcing, against the work's -users, your or third parties' legal rights to forbid circumvention of -technological measures. - - 4. Conveying Verbatim Copies. - - You may convey verbatim copies of the Program's source code as you -receive it, in any medium, provided that you conspicuously and -appropriately publish on each copy an appropriate copyright notice; -keep intact all notices stating that this License and any -non-permissive terms added in accord with section 7 apply to the code; -keep intact all notices of the absence of any warranty; and give all -recipients a copy of this License along with the Program. - - You may charge any price or no price for each copy that you convey, -and you may offer support or warranty protection for a fee. - - 5. Conveying Modified Source Versions. - - You may convey a work based on the Program, or the modifications to -produce it from the Program, in the form of source code under the -terms of section 4, provided that you also meet all of these conditions: - - a) The work must carry prominent notices stating that you modified - it, and giving a relevant date. - - b) The work must carry prominent notices stating that it is - released under this License and any conditions added under section - 7. This requirement modifies the requirement in section 4 to - "keep intact all notices". - - c) You must license the entire work, as a whole, under this - License to anyone who comes into possession of a copy. This - License will therefore apply, along with any applicable section 7 - additional terms, to the whole of the work, and all its parts, - regardless of how they are packaged. This License gives no - permission to license the work in any other way, but it does not - invalidate such permission if you have separately received it. - - d) If the work has interactive user interfaces, each must display - Appropriate Legal Notices; however, if the Program has interactive - interfaces that do not display Appropriate Legal Notices, your - work need not make them do so. - - A compilation of a covered work with other separate and independent -works, which are not by their nature extensions of the covered work, -and which are not combined with it such as to form a larger program, -in or on a volume of a storage or distribution medium, is called an -"aggregate" if the compilation and its resulting copyright are not -used to limit the access or legal rights of the compilation's users -beyond what the individual works permit. Inclusion of a covered work -in an aggregate does not cause this License to apply to the other -parts of the aggregate. - - 6. Conveying Non-Source Forms. - - You may convey a covered work in object code form under the terms -of sections 4 and 5, provided that you also convey the -machine-readable Corresponding Source under the terms of this License, -in one of these ways: - - a) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by the - Corresponding Source fixed on a durable physical medium - customarily used for software interchange. - - b) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by a - written offer, valid for at least three years and valid for as - long as you offer spare parts or customer support for that product - model, to give anyone who possesses the object code either (1) a - copy of the Corresponding Source for all the software in the - product that is covered by this License, on a durable physical - medium customarily used for software interchange, for a price no - more than your reasonable cost of physically performing this - conveying of source, or (2) access to copy the - Corresponding Source from a network server at no charge. - - c) Convey individual copies of the object code with a copy of the - written offer to provide the Corresponding Source. This - alternative is allowed only occasionally and noncommercially, and - only if you received the object code with such an offer, in accord - with subsection 6b. - - d) Convey the object code by offering access from a designated - place (gratis or for a charge), and offer equivalent access to the - Corresponding Source in the same way through the same place at no - further charge. You need not require recipients to copy the - Corresponding Source along with the object code. If the place to - copy the object code is a network server, the Corresponding Source - may be on a different server (operated by you or a third party) - that supports equivalent copying facilities, provided you maintain - clear directions next to the object code saying where to find the - Corresponding Source. Regardless of what server hosts the - Corresponding Source, you remain obligated to ensure that it is - available for as long as needed to satisfy these requirements. - - e) Convey the object code using peer-to-peer transmission, provided - you inform other peers where the object code and Corresponding - Source of the work are being offered to the general public at no - charge under subsection 6d. - - A separable portion of the object code, whose source code is excluded -from the Corresponding Source as a System Library, need not be -included in conveying the object code work. - - A "User Product" is either (1) a "consumer product", which means any -tangible personal property which is normally used for personal, family, -or household purposes, or (2) anything designed or sold for incorporation -into a dwelling. In determining whether a product is a consumer product, -doubtful cases shall be resolved in favor of coverage. For a particular -product received by a particular user, "normally used" refers to a -typical or common use of that class of product, regardless of the status -of the particular user or of the way in which the particular user -actually uses, or expects or is expected to use, the product. A product -is a consumer product regardless of whether the product has substantial -commercial, industrial or non-consumer uses, unless such uses represent -the only significant mode of use of the product. - - "Installation Information" for a User Product means any methods, -procedures, authorization keys, or other information required to install -and execute modified versions of a covered work in that User Product from -a modified version of its Corresponding Source. The information must -suffice to ensure that the continued functioning of the modified object -code is in no case prevented or interfered with solely because -modification has been made. - - If you convey an object code work under this section in, or with, or -specifically for use in, a User Product, and the conveying occurs as -part of a transaction in which the right of possession and use of the -User Product is transferred to the recipient in perpetuity or for a -fixed term (regardless of how the transaction is characterized), the -Corresponding Source conveyed under this section must be accompanied -by the Installation Information. But this requirement does not apply -if neither you nor any third party retains the ability to install -modified object code on the User Product (for example, the work has -been installed in ROM). - - The requirement to provide Installation Information does not include a -requirement to continue to provide support service, warranty, or updates -for a work that has been modified or installed by the recipient, or for -the User Product in which it has been modified or installed. Access to a -network may be denied when the modification itself materially and -adversely affects the operation of the network or violates the rules and -protocols for communication across the network. - - Corresponding Source conveyed, and Installation Information provided, -in accord with this section must be in a format that is publicly -documented (and with an implementation available to the public in -source code form), and must require no special password or key for -unpacking, reading or copying. - - 7. Additional Terms. - - "Additional permissions" are terms that supplement the terms of this -License by making exceptions from one or more of its conditions. -Additional permissions that are applicable to the entire Program shall -be treated as though they were included in this License, to the extent -that they are valid under applicable law. If additional permissions -apply only to part of the Program, that part may be used separately -under those permissions, but the entire Program remains governed by -this License without regard to the additional permissions. - - When you convey a copy of a covered work, you may at your option -remove any additional permissions from that copy, or from any part of -it. (Additional permissions may be written to require their own -removal in certain cases when you modify the work.) You may place -additional permissions on material, added by you to a covered work, -for which you have or can give appropriate copyright permission. - - Notwithstanding any other provision of this License, for material you -add to a covered work, you may (if authorized by the copyright holders of -that material) supplement the terms of this License with terms: - - a) Disclaiming warranty or limiting liability differently from the - terms of sections 15 and 16 of this License; or - - b) Requiring preservation of specified reasonable legal notices or - author attributions in that material or in the Appropriate Legal - Notices displayed by works containing it; or - - c) Prohibiting misrepresentation of the origin of that material, or - requiring that modified versions of such material be marked in - reasonable ways as different from the original version; or - - d) Limiting the use for publicity purposes of names of licensors or - authors of the material; or - - e) Declining to grant rights under trademark law for use of some - trade names, trademarks, or service marks; or - - f) Requiring indemnification of licensors and authors of that - material by anyone who conveys the material (or modified versions of - it) with contractual assumptions of liability to the recipient, for - any liability that these contractual assumptions directly impose on - those licensors and authors. - - All other non-permissive additional terms are considered "further -restrictions" within the meaning of section 10. If the Program as you -received it, or any part of it, contains a notice stating that it is -governed by this License along with a term that is a further -restriction, you may remove that term. If a license document contains -a further restriction but permits relicensing or conveying under this -License, you may add to a covered work material governed by the terms -of that license document, provided that the further restriction does -not survive such relicensing or conveying. - - If you add terms to a covered work in accord with this section, you -must place, in the relevant source files, a statement of the -additional terms that apply to those files, or a notice indicating -where to find the applicable terms. - - Additional terms, permissive or non-permissive, may be stated in the -form of a separately written license, or stated as exceptions; -the above requirements apply either way. - - 8. Termination. - - You may not propagate or modify a covered work except as expressly -provided under this License. Any attempt otherwise to propagate or -modify it is void, and will automatically terminate your rights under -this License (including any patent licenses granted under the third -paragraph of section 11). - - However, if you cease all violation of this License, then your -license from a particular copyright holder is reinstated (a) -provisionally, unless and until the copyright holder explicitly and -finally terminates your license, and (b) permanently, if the copyright -holder fails to notify you of the violation by some reasonable means -prior to 60 days after the cessation. - - Moreover, your license from a particular copyright holder is -reinstated permanently if the copyright holder notifies you of the -violation by some reasonable means, this is the first time you have -received notice of violation of this License (for any work) from that -copyright holder, and you cure the violation prior to 30 days after -your receipt of the notice. - - Termination of your rights under this section does not terminate the -licenses of parties who have received copies or rights from you under -this License. If your rights have been terminated and not permanently -reinstated, you do not qualify to receive new licenses for the same -material under section 10. - - 9. Acceptance Not Required for Having Copies. - - You are not required to accept this License in order to receive or -run a copy of the Program. Ancillary propagation of a covered work -occurring solely as a consequence of using peer-to-peer transmission -to receive a copy likewise does not require acceptance. However, -nothing other than this License grants you permission to propagate or -modify any covered work. These actions infringe copyright if you do -not accept this License. Therefore, by modifying or propagating a -covered work, you indicate your acceptance of this License to do so. - - 10. Automatic Licensing of Downstream Recipients. - - Each time you convey a covered work, the recipient automatically -receives a license from the original licensors, to run, modify and -propagate that work, subject to this License. You are not responsible -for enforcing compliance by third parties with this License. - - An "entity transaction" is a transaction transferring control of an -organization, or substantially all assets of one, or subdividing an -organization, or merging organizations. If propagation of a covered -work results from an entity transaction, each party to that -transaction who receives a copy of the work also receives whatever -licenses to the work the party's predecessor in interest had or could -give under the previous paragraph, plus a right to possession of the -Corresponding Source of the work from the predecessor in interest, if -the predecessor has it or can get it with reasonable efforts. - - You may not impose any further restrictions on the exercise of the -rights granted or affirmed under this License. For example, you may -not impose a license fee, royalty, or other charge for exercise of -rights granted under this License, and you may not initiate litigation -(including a cross-claim or counterclaim in a lawsuit) alleging that -any patent claim is infringed by making, using, selling, offering for -sale, or importing the Program or any portion of it. - - 11. Patents. - - A "contributor" is a copyright holder who authorizes use under this -License of the Program or a work on which the Program is based. The -work thus licensed is called the contributor's "contributor version". - - A contributor's "essential patent claims" are all patent claims -owned or controlled by the contributor, whether already acquired or -hereafter acquired, that would be infringed by some manner, permitted -by this License, of making, using, or selling its contributor version, -but do not include claims that would be infringed only as a -consequence of further modification of the contributor version. For -purposes of this definition, "control" includes the right to grant -patent sublicenses in a manner consistent with the requirements of -this License. - - Each contributor grants you a non-exclusive, worldwide, royalty-free -patent license under the contributor's essential patent claims, to -make, use, sell, offer for sale, import and otherwise run, modify and -propagate the contents of its contributor version. - - In the following three paragraphs, a "patent license" is any express -agreement or commitment, however denominated, not to enforce a patent -(such as an express permission to practice a patent or covenant not to -sue for patent infringement). To "grant" such a patent license to a -party means to make such an agreement or commitment not to enforce a -patent against the party. - - If you convey a covered work, knowingly relying on a patent license, -and the Corresponding Source of the work is not available for anyone -to copy, free of charge and under the terms of this License, through a -publicly available network server or other readily accessible means, -then you must either (1) cause the Corresponding Source to be so -available, or (2) arrange to deprive yourself of the benefit of the -patent license for this particular work, or (3) arrange, in a manner -consistent with the requirements of this License, to extend the patent -license to downstream recipients. "Knowingly relying" means you have -actual knowledge that, but for the patent license, your conveying the -covered work in a country, or your recipient's use of the covered work -in a country, would infringe one or more identifiable patents in that -country that you have reason to believe are valid. - - If, pursuant to or in connection with a single transaction or -arrangement, you convey, or propagate by procuring conveyance of, a -covered work, and grant a patent license to some of the parties -receiving the covered work authorizing them to use, propagate, modify -or convey a specific copy of the covered work, then the patent license -you grant is automatically extended to all recipients of the covered -work and works based on it. - - A patent license is "discriminatory" if it does not include within -the scope of its coverage, prohibits the exercise of, or is -conditioned on the non-exercise of one or more of the rights that are -specifically granted under this License. You may not convey a covered -work if you are a party to an arrangement with a third party that is -in the business of distributing software, under which you make payment -to the third party based on the extent of your activity of conveying -the work, and under which the third party grants, to any of the -parties who would receive the covered work from you, a discriminatory -patent license (a) in connection with copies of the covered work -conveyed by you (or copies made from those copies), or (b) primarily -for and in connection with specific products or compilations that -contain the covered work, unless you entered into that arrangement, -or that patent license was granted, prior to 28 March 2007. - - Nothing in this License shall be construed as excluding or limiting -any implied license or other defenses to infringement that may -otherwise be available to you under applicable patent law. - - 12. No Surrender of Others' Freedom. - - If conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot convey a -covered work so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you may -not convey it at all. For example, if you agree to terms that obligate you -to collect a royalty for further conveying from those to whom you convey -the Program, the only way you could satisfy both those terms and this -License would be to refrain entirely from conveying the Program. - - 13. Use with the GNU Affero General Public License. - - Notwithstanding any other provision of this License, you have -permission to link or combine any covered work with a work licensed -under version 3 of the GNU Affero General Public License into a single -combined work, and to convey the resulting work. The terms of this -License will continue to apply to the part which is the covered work, -but the special requirements of the GNU Affero General Public License, -section 13, concerning interaction through a network will apply to the -combination as such. - - 14. Revised Versions of this License. - - The Free Software Foundation may publish revised and/or new versions of -the GNU General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - - Each version is given a distinguishing version number. If the -Program specifies that a certain numbered version of the GNU General -Public License "or any later version" applies to it, you have the -option of following the terms and conditions either of that numbered -version or of any later version published by the Free Software -Foundation. If the Program does not specify a version number of the -GNU General Public License, you may choose any version ever published -by the Free Software Foundation. - - If the Program specifies that a proxy can decide which future -versions of the GNU General Public License can be used, that proxy's -public statement of acceptance of a version permanently authorizes you -to choose that version for the Program. - - Later license versions may give you additional or different -permissions. However, no additional obligations are imposed on any -author or copyright holder as a result of your choosing to follow a -later version. - - 15. Disclaimer of Warranty. - - THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY -APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT -HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY -OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, -THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM -IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF -ALL NECESSARY SERVICING, REPAIR OR CORRECTION. - - 16. Limitation of Liability. - - IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS -THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY -GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE -USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF -DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD -PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), -EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF -SUCH DAMAGES. - - 17. Interpretation of Sections 15 and 16. - - If the disclaimer of warranty and limitation of liability provided -above cannot be given local legal effect according to their terms, -reviewing courts shall apply local law that most closely approximates -an absolute waiver of all civil liability in connection with the -Program, unless a warranty or assumption of liability accompanies a -copy of the Program in return for a fee. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -state the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - - Copyright (C) - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . - -Also add information on how to contact you by electronic and paper mail. +Copyright Christian Hoener zu Siederdissen 2011-2019 - If the program does terminal interaction, make it output a short -notice like this when it starts in an interactive mode: +All rights reserved. - Copyright (C) - This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. - This is free software, and you are welcome to redistribute it - under certain conditions; type `show c' for details. +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: -The hypothetical commands `show w' and `show c' should show the appropriate -parts of the General Public License. Of course, your program's commands -might be different; for a GUI interface, you would use an "about box". + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. - You should also get your employer (if you work as a programmer) or school, -if any, to sign a "copyright disclaimer" for the program, if necessary. -For more information on this, and how to apply and follow the GNU GPL, see -. + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. - The GNU General Public License does not permit incorporating your program -into proprietary programs. If your program is a subroutine library, you -may consider it more useful to permit linking proprietary applications with -the library. If this is what you want to do, use the GNU Lesser General -Public License instead of this License. But first, please read -. + * Neither the name of Christian Hoener zu Siederdissen nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,27 @@ +![github action: master](https://github.com/choener/BiobaseXNA/actions/workflows/ci.yml/badge.svg?branch=master) +![github action: hackage](https://github.com/choener/SciBaseTypes/actions/workflows/hackage.yml/badge.svg) + +# BiobaseXNA + +Efficient encoding of (short) biological sequences. This package ist designed +to deal with *in-memory* snippets of DNA, RNA, and amino acids. The encoding is +geared toward time-efficiency, not necessarily space efficiency (we use Int's +for encoding characters, not the smallest type possible). + +Additional modules provide conversion capabilities between different types of +characters according to biological laws, and some biochemical constraint +information. The latter includes canonical and non-canonical pairing +information for RNA. + +Actual energy parameters for pairings are provided by other packages, for +example BiobaseTurner for the loop energy model with measured parameters. + + + +#### Contact + +Christian Hoener zu Siederdissen +Leipzig University, Leipzig, Germany +choener@bioinf.uni-leipzig.de +http://www.bioinf.uni-leipzig.de/~choener/ + diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,92 @@ +0.11.0.1 +-------- + +- version bump + +0.11.0.0 +-------- + +0.10.0.0 +-------- + +- redesigned Biobase.Secondary.Basepair + +0.9.3.1 +------- + +- removed upper version bounds, bumped dependent versions +- added d1Distance for fast basepair distance calculations + +0.9.3.0 +------- + +- bigger version bump together with multiple ghc compatibility + +0.9.2.1 +------- + +- stack.yaml, some version bumping + +0.9.2.0 +------- + +- fixed overlapping instances in AA + +0.9.1.0 +------- + +- introduction of pattern synonyms for letters in molecular alphabets +- travis-ci integration + +0.9.0.0 +------- + +- major cleanup of the XNA library: explicit encoding of RNA,DNA, or XNA (XNA + contains both U and T) +- IUPAC notation (degenerate nucleotides) has an efficient encoding as well +- translation between XNA and protein (Biobase.Primary.Trans) +- import Biobase.Primary to get everything under Primary.* +- import Biobase.Secondary to get everything under Secondary.* +- SubOptDistance now extends all structure lines of RNAsubopt with a third + field, the distance +- Diagrams provide methods to validate folding and cofolding secondary + structure strings. +- serialization capabilities for 'Letter's + +0.8.3.0 +------- + +- bugfix version: use vector-th-unbox to generate unboxed vector instances + +0.8.2.0 +------- + +- dotBracket -> unsafeDotBracket +- new 'dotBracket' function works in the error monad + +0.8.1.1 +------- + +- added T/U conversion functions + +0.8.1.0 +------- + +- Biobase.Primary.IUPAC for degenerate base symbol conversion + +0.8.0.0 +------- + +- Biobase.Codon -> Biobase.AAseq +- and efficient encoding of AAseqs + +0.7 +--- + +- updated to PrimitiveArray >= 0.5 +- added Codon table + +0.6.2.0 +------- + +- Updated to PrimitiveArray >= 0.2.0.0 diff --git a/sources/codontable b/sources/codontable new file mode 100644 --- /dev/null +++ b/sources/codontable @@ -0,0 +1,64 @@ +AAA K +AAC N +AAG K +AAT N +ACA T +ACC T +ACG T +ACT T +AGA R +AGC S +AGG R +AGT S +ATA I +ATC I +ATG M +ATT I +CAA Q +CAC H +CAG Q +CAT H +CCA P +CCC P +CCG P +CCT P +CGA R +CGC R +CGG R +CGT R +CTA L +CTC L +CTG L +CTT L +GAA E +GAC D +GAG E +GAT D +GCA A +GCC A +GCG A +GCT A +GGA G +GGC G +GGG G +GGT G +GTA V +GTC V +GTG V +GTT V +TAA * +TAC Y +TAG * +TAT Y +TCA S +TCC S +TCG S +TCT S +TGA * +TGC C +TGG W +TGT C +TTA L +TTC F +TTG L +TTT F diff --git a/sources/iupac-nucleotides b/sources/iupac-nucleotides new file mode 100644 --- /dev/null +++ b/sources/iupac-nucleotides @@ -0,0 +1,16 @@ +A A +C C +G G +T T +U U +W AT +S CG +M AC +K GT +R AG +Y CT +B CGT +D AGT +H ACT +V ACG +N ACGT diff --git a/src/SubOptDistance.hs b/src/SubOptDistance.hs new file mode 100644 --- /dev/null +++ b/src/SubOptDistance.hs @@ -0,0 +1,45 @@ + +-- | A very simple program that calculates base pair distances given +-- a secondary structure as argument and a Vienna RNAsubopt output via +-- stdin. + +module Main where + +import Control.Arrow +import Data.Char (isSpace) +import Data.Either (either) +import Data.Ord +import System.Console.CmdArgs +import Text.Printf + +import Biobase.Secondary +import Biobase.Secondary.Diagrams + + + +data Options = Options + { structure :: String + , aq :: Bool + , qa :: Bool + } deriving (Show,Data,Typeable) + +options = Options + { structure = "" &= args + , aq = True &= help "perform target \\\\ query calculation" + , qa = True &= help "perform query \\\\ target calculation" + } &= help helpLines + +helpLines = unlines + [ "This program reads RNAsubopt output on STDIN. Provide an" + , "additional structure as the argument line. The result will be the" + , "sub-optimal line with lowest base pair distance to the query line." + ] + +main = do + Options{..} <- cmdArgs options + (sqn:xs') <- fmap lines $ getContents + putStrLn sqn + let s' = fst . break isSpace $ structure + let xs = map (viennaStringDistance aq qa s') xs' + mapM_ (\(s,d) -> printf "%s %6d\n" s d) xs + diff --git a/tests/properties.hs b/tests/properties.hs new file mode 100644 --- /dev/null +++ b/tests/properties.hs @@ -0,0 +1,116 @@ + +module Main where + +import Control.Monad (join) +import Data.Bits +import Data.Function (on) +import Data.Int (Int16(..)) +import Data.List (groupBy,sort,permutations,nub,(\\)) +import Data.Maybe (isJust) +import Data.Word (Word) +import Debug.Trace +import qualified Data.Vector.Unboxed as VU +import Test.QuickCheck hiding ((.&.)) +import Test.Tasty +import Test.Tasty.QuickCheck (testProperty) +import Test.Tasty.TH + +import Biobase.Secondary.Diagrams +import Biobase.Secondary.Basepair (PairIdx) + + + +newtype ArbitrarySSTree = ASST (SSTree PairIdx ()) + deriving (Show) + +instance Arbitrary ArbitrarySSTree where + arbitrary = ASST <$> sized arbitrarySSTree + where + arbitrarySSTree m = do + Positive c <- arbitrary + cs <- go 0 (c*m) + let k = if null cs then 0 else 1 + maximum [ z | SSTree (_,z) _ _ <- cs ] + return $ SSExtern k () cs + go i j = do + Positive c <- arbitrary + Positive d <- arbitrary + if i+c+d >= j + then return [] + else do + cs <- go (i+c+1) (i+c+d) + let h = SSTree (i+c,i+c+d) () cs + ts <- go (i+c+d+1) j + return $ h:ts +{- + shrink (ASST (SSExtern k () cs)) + | null cs = [] + | otherwise = [ ASST +-} + +collectPairs (ASST (SSExtern k _ zs)) = (k, sort $ go zs) + where go [] = [] + go (SSTree (i,j) _ cs : ss) = (i,j) : go cs ++ go ss + +bld :: Int -> [PairIdx] -> D1Secondary +bld = curry mkD1S + +prop_d1Distance a@(ASST _) b@(ASST _) = d1Distance x y == k + where x = mkD1S (lx', x') + y = mkD1S (ly', y') + (lx', x') = collectPairs a + (ly', y') = collectPairs b + k = length $ (x' \\ y') ++ (y' \\ x') + +---- | Check if both the memoized version and the population enumeration +---- produce the same multisets, but maybe in different order. +---- +---- prop> \(n :: Int16) -> let b = popCount n in memoSorted b == enumSorted b +---- +-- +--prop_PopCountSet (NonZero (n' :: Int16)) = memo == enum +-- where b = popCount n +-- memo = memoSorted b +-- enum = enumSorted b +-- n = n' `mod` 12 +-- +--memoSorted, enumSorted :: Int -> [[Int]] +-- +--memoSorted b = map sort . groupBy ((==) `on` popCount) $ VU.toList $ popCntMemoInt b +--enumSorted b = map sort $ [0] : [ roll (popPermutation b) (Just $ 2^k-1) | k <- [1..b] ] +-- where roll f (Just k) = k : roll f (f k) +-- roll _ Nothing = [] +-- +--prop_lsb_Int (x :: Int) = lsbZ x == maybe (-1) id (maybeLsb x) +-- +--prop_lsb_Word (x :: Word) = lsbZ x == maybe (-1) id (maybeLsb x) +-- +--prop_OneBits_Int (x :: Int) = popCount x == length abl && and [ testBit x k | k <- abl ] +-- where abl = activeBitsL x +-- +---- Tests if we actually generate all permutations. +-- +--prop_allPermutations (a :: Int , b :: Int) = and $ zipWith cmp (sort qs) (sort $ nub ps) +-- where nbs = min a' b' -- number of 1 bits in set +-- sts = max a' b' -- set size +-- a' = a `mod` 8 -- finiteBitSize a +-- b' = b `mod` 8 -- finiteBitSize b +-- ps = permutations $ replicate (sts - nbs) False ++ replicate nbs True +-- qs = go (Just $ 2 ^ nbs - 1) +-- go :: Maybe Int -> [Int] +-- go Nothing = [] +-- go (Just k) = k : go (popPermutation sts k) +-- cmp k as = and [ if a then testBit k c else (not $ testBit k c) | (a,c) <- zip (reverse as) [0 .. ] ] +-- +---- TODO popComplement +-- +--prop_popShiftL_popShiftR (a::Word,b::Word) = s == l +-- where m = a .|. b +-- s = a .&. b +-- l = popShiftL m r +-- r = popShiftR m s + + + +main :: IO () +main = $(defaultMainGenerator) +