diff --git a/Biobase/AAseq.hs b/Biobase/AAseq.hs deleted file mode 100644 --- a/Biobase/AAseq.hs +++ /dev/null @@ -1,254 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE ViewPatterns #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE TypeFamilies #-} - --- | This module has the translation tables for the genetic code. - -module Biobase.AAseq where - -import Control.Arrow ((***)) -import Data.Ix (Ix(..)) -import Data.Primitive.Types -import Data.Tuple (swap) -import Data.Vector.Unboxed.Deriving -import GHC.Base (remInt,quotInt) -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 Data.Array.Repa.ExtShape -import Data.Array.Repa.Index -import Data.Array.Repa.Shape - -import Biobase.Primary - - - --- | The amino acid newtype. - -newtype AA = AA { unAA :: Int } - deriving (Eq,Ord,Ix) - - - --- * Creating functions and aa data. - -(aStop:aA:aB:aC:aD:aE:aF:aG:aH:aI:aK:aL:aM:aN:aP:aQ:aR:aS:aT:aV:aW:aX:aY:aZ:aUndefined:_) = map AA [0..] - -aaRange = [aStop .. pred aUndefined] - --- | Translate 'Char' amino acid representation into efficient 'AA' newtype. - -toAA :: Char -> AA -toAA ((`lookup` charAAList) -> Just aa) = aa -toAA c = error $ "unknown AA: " ++ show c - --- | 'Char' representation of an 'AA'. - -fromAA :: AA -> Char -fromAA ((`lookup` aACharList) -> Just c) = c -fromAA (AA aa) = error $ "unknown AA: " ++ (show aa) - --- | Create amino acid sequences from different data sources. - -class MkAAseq x where - mkAAseq :: x -> VU.Vector AA - -type AAseq = VU.Vector AA - --- | Using the codon table, create an 'AAseq' from the 'Primary' sequence. - -primaryToAAseq :: Primary -> AAseq -primaryToAAseq = mkAAseq . go where - go (VU.length -> 0) = [] - go (VU.splitAt 3 -> (VU.toList -> hs,ts)) = case M.lookup hs nucCodonTable of - Just aa -> aa : go ts - _ -> error $ "primaryToAAseq: " ++ show (hs,ts) - - - --- * lookup tables - -charAAList = - [ ('/',aStop) - , ('A',aA) - , ('B',aB) - , ('C',aC) - , ('D',aD) - , ('E',aE) - , ('F',aF) - , ('G',aG) - , ('H',aH) - , ('I',aI) - , ('K',aK) - , ('L',aL) - , ('M',aM) - , ('N',aN) - , ('P',aP) - , ('Q',aQ) - , ('R',aR) - , ('S',aS) - , ('T',aT) - , ('V',aV) - , ('W',aW) - , ('X',aX) - , ('Y',aY) - , ('Z',aZ) - ] - -aACharList = map swap charAAList - -codonTable = M.fromList - [ ("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') - ] - -nucCodonTable = M.fromList . map (map mkNuc *** toAA) . M.assocs $ codonTable - - - --- * instances - -instance Show AA where - show n = [fromAA n] - -instance Read AA where - readsPrec p [] = [] - readsPrec p (x:xs) - | x==' ' = readsPrec p xs - | Just aa <- x `lookup` charAAList = [(aa,xs)] - | otherwise = [] - -instance (Shape sh,Show sh) => Shape (sh :. AA) where - rank (sh:._) = rank sh + 1 - zeroDim = zeroDim:.AA 0 - unitDim = unitDim:.AA 1 -- TODO does this one make sense? - intersectDim (sh1:.n1) (sh2:.n2) = intersectDim sh1 sh2 :. min n1 n2 - addDim (sh1:.AA n1) (sh2:.AA n2) = addDim sh1 sh2 :. AA (n1+n2) -- TODO will not necessarily yield a valid Nuc - size (sh1:.AA n) = size sh1 * n - sizeIsValid (sh1:.AA n) = sizeIsValid (sh1:.n) - toIndex (sh1:.AA sh2) (sh1':.AA sh2') = toIndex (sh1:.sh2) (sh1':.sh2') - fromIndex (ds:.AA d) n = fromIndex ds (n `quotInt` d) :. AA 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 :. AA 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 :. AA) where - subDim (sh1:.AA n1) (sh2:.AA n2) = subDim sh1 sh2 :. AA (n1-n2) - rangeList (sh1:.AA n1) (sh2:.AA n2) = [ sh:.AA n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2)]] - -instance Enum AA where - toEnum = AA - fromEnum = unAA - -instance MkAAseq [Char] where - mkAAseq = VU.fromList . map toAA - -instance MkAAseq [AA] where - mkAAseq = VU.fromList - -instance MkAAseq (VU.Vector Char) where - mkAAseq = VU.map toAA - -instance MkAAseq BS.ByteString where - mkAAseq = VU.fromList . map toAA . BS.unpack - -instance MkAAseq BSL.ByteString where - mkAAseq = VU.fromList . map toAA . BSL.unpack - -instance MkAAseq T.Text where - mkAAseq = VU.fromList . map toAA . T.unpack - -derivingUnbox "AA" - [t| AA -> Int |] [| unAA |] [| AA |] - diff --git a/Biobase/Primary.hs b/Biobase/Primary.hs --- a/Biobase/Primary.hs +++ b/Biobase/Primary.hs @@ -1,205 +1,21 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE OverlappingInstances #-} -{-# LANGUAGE PackageImports #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE TypeSynonymInstances #-} --- | 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 Data.Char (toUpper) -import Data.Ix (Ix(..)) -import Data.Primitive.Types -import Data.Tuple (swap) -import Data.Vector.Unboxed.Deriving -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.Array.Repa.ExtShape -import Data.Array.Repa.Index -import Data.Array.Repa.Shape - -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 = VU.Vector Nuc - - - - --- * 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 : nU : nIMI : nUndefined : _) = map Nuc [0 .. ] - -acgt = [nA,nC,nG,nT] -acgu = [nA,nC,nG,nU] -cgau = [nC,nG,nA,nU] -nacgt = nN:acgt -nacgu = nN:acgu - --- | 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 - --- ** Methods to convert between DNA and RNA +-- | -- --- TODO add all the rev-comp stuff and whatnot - --- | @T@ to @U@ - -convT2U x - | x == nT = nU - | otherwise = x - --- | @U@ to @T@ - -convU2T x - | x == nU = nT - | otherwise = x - - - --- * 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 = [] - -derivingUnbox "Nuc" - [t| Nuc -> Int |] [| unNuc |] [| 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)]] - --- | Enum - -instance Enum Nuc where - toEnum = Nuc - fromEnum = unNuc - --- ** Instances for 'MkPrimary' - -instance MkPrimary String where - mkPrimary = VU.fromList . map mkNuc - -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 + ) where -instance MkPrimary [Nuc] where - mkPrimary = VU.fromList +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 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,152 @@ + +-- | This module has the translation tables for the genetic code. We do +-- have a symbol 'Undef' for undefined amino acids (say because of @N@s in +-- the nucleotide code). + +module Biobase.Primary.AA where + +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 GHC.Base (remInt,quotInt) +import GHC.Generics (Generic) +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 Data.Bijection.Map as B + +import Biobase.Primary.Letter + + + +-- | Amino acid phantom type. + +data AA + +pattern Stop = Letter 0 :: Letter AA +pattern A = Letter 1 :: Letter AA +pattern B = Letter 2 :: Letter AA +pattern C = Letter 3 :: Letter AA +pattern D = Letter 4 :: Letter AA +pattern E = Letter 5 :: Letter AA +pattern F = Letter 6 :: Letter AA +pattern G = Letter 7 :: Letter AA +pattern H = Letter 8 :: Letter AA +pattern I = Letter 9 :: Letter AA +pattern K = Letter 10 :: Letter AA +pattern L = Letter 11 :: Letter AA +pattern M = Letter 12 :: Letter AA +pattern N = Letter 13 :: Letter AA +pattern P = Letter 14 :: Letter AA +pattern Q = Letter 15 :: Letter AA +pattern R = Letter 16 :: Letter AA +pattern S = Letter 17 :: Letter AA +pattern T = Letter 18 :: Letter AA +pattern V = Letter 19 :: Letter AA +pattern W = Letter 20 :: Letter AA +pattern X = Letter 21 :: Letter AA +pattern Y = Letter 22 :: Letter AA +pattern Z = Letter 23 :: Letter AA +pattern Undef = Letter 24 :: Letter AA + + +-- * Creating functions and aa data. + +aa :: Int -> Letter AA +aa = Letter + +aaRange = [Stop .. pred Undef] + +-- | Translate 'Char' amino acid representation into efficient 'AA' newtype. + +charAA :: Char -> Letter AA +charAA = B.findWithDefaultL Undef charBaa +{-# INLINE charAA #-} + +-- | 'Char' representation of an 'AA'. + +aaChar :: Letter AA -> Char +aaChar = B.findWithDefaultR '?' charBaa +{-# INLINE aaChar #-} + +-- * lookup tables + +charBaa :: B.Bimap Char (Letter AA) +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) + , ('?',Undef) + ] +{-# NOINLINE charBaa #-} + + + +-- * instances + +instance Show (Letter AA) where + show n = [aaChar n] + +instance Read (Letter AA) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | aa <- charAA x = [(aa,xs)] + | otherwise = [] + +instance Enum (Letter AA) 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<=(unLetter Undef) = Letter k + toEnum k = error $ "toEnum/Letter RNA " ++ show k + fromEnum (Letter k) = k + +instance MkPrimary [Char] AA where + primary = VU.fromList . map charAA + +instance MkPrimary [Letter AA] AA where + primary = VU.fromList + +instance MkPrimary (VU.Vector Char) AA where + primary = VU.map charAA + +instance MkPrimary BS.ByteString AA where + primary = VU.fromList . map charAA . BS.unpack + +instance MkPrimary BSL.ByteString AA where + primary = VU.fromList . map charAA . BSL.unpack + +instance MkPrimary T.Text AA where + primary = VU.fromList . map charAA . T.unpack + diff --git a/Biobase/Primary/Bounds.hs b/Biobase/Primary/Bounds.hs deleted file mode 100644 --- a/Biobase/Primary/Bounds.hs +++ /dev/null @@ -1,52 +0,0 @@ - --- | A special class of bounds for RNA/pair encodings that are used to index --- into tables. We typically encode more in the alphabets than we want to use --- to index, so in order to keep things simple, we have specialized bounds. - -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 --- Biobase.RNA.ViennaPair for instances. - -class Bounded a => Bounds a where - minNormal :: a - maxNormal :: a - minExtended :: a - maxExtended :: a - - - --- * Instances for tuples of size 2-6 - -instance (Bounds a, Bounds b) => Bounds (a,b) where - minNormal = (minNormal, minNormal) - maxNormal = (maxNormal, maxNormal) - minExtended = (minExtended, minExtended) - maxExtended = (maxExtended, maxExtended) - -instance (Bounds a, Bounds b, Bounds c) => Bounds (a,b,c) where - minNormal = (minNormal, minNormal, minNormal) - maxNormal = (maxNormal, maxNormal, maxNormal) - minExtended = (minExtended, minExtended, minExtended) - maxExtended = (maxExtended, maxExtended, maxExtended) - -instance (Bounds a, Bounds b, Bounds c, Bounds d) => Bounds (a,b,c,d) where - minNormal = (minNormal, minNormal, minNormal, minNormal) - maxNormal = (maxNormal, maxNormal, maxNormal, maxNormal) - minExtended = (minExtended, minExtended, minExtended, minExtended) - maxExtended = (maxExtended, maxExtended, maxExtended, maxExtended) - -instance (Bounds a, Bounds b, Bounds c, Bounds d, Bounds e) => Bounds (a,b,c,d,e) where - minNormal = (minNormal, minNormal, minNormal, minNormal, minNormal) - maxNormal = (maxNormal, maxNormal, maxNormal, maxNormal, maxNormal) - minExtended = (minExtended, minExtended, minExtended, minExtended, minExtended) - maxExtended = (maxExtended, maxExtended, maxExtended, maxExtended, maxExtended) - -instance (Bounds a, Bounds b, Bounds c, Bounds d, Bounds e, Bounds f) => Bounds (a,b,c,d,e,f) where - minNormal = (minNormal, minNormal, minNormal, minNormal, minNormal, minNormal) - 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,17 +1,11 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE TypeFamilies #-} --- | Fast hash functions for 'Primary' sequences. A hash is just an 'Int', so --- use these only for short sequences. --- --- TODO replace with standard hashing functions used by Haskell libs? +-- | 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.Vector.Unboxed.Deriving @@ -19,34 +13,34 @@ import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU -import Biobase.Primary +import Biobase.Primary.Letter -newtype HashedPrimary = HashedPrimary { unHashedPrimary :: Int } +-- | The hash of a primary sequence. + +newtype HashedPrimary t = HashedPrimary { unHashedPrimary :: Int } deriving (Eq,Ord,Ix,Read,Show,Enum,Bounded) derivingUnbox "HashedPrimary" - [t| HashedPrimary -> Int |] [| unHashedPrimary |] [| HashedPrimary |] - + [t| forall a . HashedPrimary a -> 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 +mkHashedPrimary :: forall t . (VU.Unbox (Letter t), Bounded (Letter t), Enum (Letter t)) => Primary t -> HashedPrimary t +mkHashedPrimary = HashedPrimary . fst . VU.foldl' f (0, 1) where + f (z, c) n = (z + c * (fromEnum n +1), c * (fromEnum (maxBound :: Letter t) + 1)) {-# INLINE mkHashedPrimary #-} + +-- | Turn a hash back into a sequence. Will fail if the resulting sequence +-- has more than 100 elements. + +hash2primary :: forall t . (VU.Unbox (Letter t), Bounded (Letter t), Enum (Letter t)) => HashedPrimary t -> Primary t +hash2primary (HashedPrimary h) = VU.unfoldrN l f h where + m = fromEnum (maxBound :: Letter t) +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 --- a/Biobase/Primary/IUPAC.hs +++ b/Biobase/Primary/IUPAC.hs @@ -1,47 +1,150 @@ -{-# LANGUAGE TemplateHaskell #-} -- | Degenerate base symbol representation. We use the same conventions as in -- <> which ignores --- @U@racil, except if it stands alone. Therefore, any RNA sequence should be --- converted to DNA (and back afterwards). --- --- NOTE that the generic 'Char' instance is not optimized for speed. +-- @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 this should be easier once we have instances for RNA,DNA, etc +-- TODO Shall we handle 'Complement' for degenerates? module Biobase.Primary.IUPAC where -import Data.ByteString.Char8 (ByteString,unpack) -import Data.FileEmbed (embedFile) -import Data.List (nub,sort) -import Data.Tuple (swap) +import Control.Arrow ((***)) +import Data.ByteString.Char8 (ByteString,unpack) +import Data.Char (toUpper) +import Data.FileEmbed (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.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 +pattern C = Letter 1 :: Letter DEG +pattern G = Letter 2 :: Letter DEG +pattern T = Letter 3 :: Letter DEG +pattern U = Letter 4 :: Letter DEG +pattern W = Letter 5 :: Letter DEG +pattern S = Letter 6 :: Letter DEG +pattern M = Letter 7 :: Letter DEG +pattern K = Letter 8 :: Letter DEG +pattern R = Letter 9 :: Letter DEG +pattern Y = Letter 10 :: Letter DEG +pattern B = Letter 11 :: Letter DEG +pattern D = Letter 12 :: Letter DEG +pattern H = Letter 13 :: Letter DEG +pattern V = Letter 14 :: Letter DEG +pattern N = Letter 15 :: Letter DEG + +instance Bounded (Letter DEG) where + minBound = A + maxBound = N + +instance Enum (Letter DEG) 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) where + show c = [degChar c] + +degSeq :: MkPrimary n DEG => n -> Primary DEG +degSeq = primary + +instance MkPrimary (VU.Vector Char) DEG where + primary = VU.map charDEG + +instance IsString [Letter DEG] where + fromString = map charDEG + + + +-- * Conversions + class Degenerate x where - fromSymbol :: x -> [x] - toSymbol :: [x] -> Maybe x + fromDegenerate :: Char -> [x] + toDegenerate :: [x] -> Maybe Char instance Degenerate Char where - fromSymbol = maybe [] id . flip lookup iupacList - toSymbol = flip lookup (map swap iupacList) . nub . sort + fromDegenerate = maybe [] id . flip lookup iupacXDNAchars + toDegenerate = flip lookup (map swap iupacXDNAchars) . nub . sort --- instance Degenerate RNA where --- --- instance Degenerate DNA where --- --- instance Degenerate XNA where -- if we want a combined alphabet +instance Degenerate (Letter RNA) 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) where + fromDegenerate 'U' = [] + fromDegenerate x = map charDNA $ fromDegenerate x + toDegenerate = toDegenerate . map dnaChar +instance Degenerate (Letter XNA) where + fromDegenerate = map charXNA . fromDegenerate + toDegenerate = toDegenerate . map xnaChar --- ** Raw embeddings --- | list of characters -iupacList :: [(Char,String)] -iupacList = map (go . words) . lines . unpack $ iupacNucleotides where +-- * 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 iupacList #-} +{-# NOINLINE iupacXDNAchars #-} -- | Raw iupac data, embedded into the library. 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,158 @@ + +-- | A newtype with an attached phenotype 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 Data.Aeson +import Data.Binary +import Data.Hashable (Hashable) +import Data.Ix (Ix(..)) +import Data.Serialize (Serialize(..)) +import Data.String (IsString(..)) +import Data.Vector.Fusion.Stream.Monadic (map,flatten,Step(..)) +import Data.Vector.Fusion.Stream.Size (Size (Unknown)) +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) + + + +-- | A 'Letter' together with its phantom type @t@ encodes bio-sequences. + +newtype Letter t = Letter { unLetter :: Int } + deriving (Eq,Ord,Generic,Ix) + +instance Binary (Letter t) +instance Serialize (Letter t) +instance FromJSON (Letter t) +instance ToJSON (Letter t) + +type Primary t = VU.Vector (Letter t) + +-- | Conversion from a large number of sequence-like inputs to primary +-- sequences. + +class MkPrimary n t where + primary :: n -> Primary t + +instance (MkPrimary (VU.Vector Char) t) => MkPrimary String t where + primary = primary . VU.fromList + +instance MkPrimary (VU.Vector Char) t => MkPrimary T.Text t where + primary = primary . VU.fromList . T.unpack + +instance MkPrimary (VU.Vector Char) t => MkPrimary TL.Text t where + primary = primary . VU.fromList . TL.unpack + +instance MkPrimary (VU.Vector Char) t => MkPrimary BS.ByteString t where + primary = primary . VU.fromList . BS.unpack + +instance MkPrimary (VU.Vector Char) t => MkPrimary BSL.ByteString t where + primary = primary . VU.fromList . BSL.unpack + +instance (VU.Unbox (Letter t), IsString [Letter t]) => IsString (VU.Vector (Letter t)) where + fromString = VU.fromList . fromString + + + +-- *** Instances for 'Letter'. + +derivingUnbox "Letter" + [t| forall a . Letter a -> Int |] [| unLetter |] [| Letter |] + +instance Hashable (Letter t) + +instance Index (Letter l) where + linearIndex _ _ (Letter i) = i + {-# Inline linearIndex #-} + smallestLinearIndex _ = error "still needed?" + {-# Inline smallestLinearIndex #-} + largestLinearIndex (Letter h) = h + {-# Inline largestLinearIndex #-} + size _ (Letter h) = h+1 + {-# Inline size #-} + inBounds (Letter l) (Letter h) (Letter i) = l <= i && i <= h + {-# Inline inBounds #-} + +instance IndexStream z => IndexStream (z:.Letter l) where + streamUp (ls:.Letter l) (hs:.Letter h) = flatten mk step Unknown $ streamUp ls hs + where mk z = return (z,l) + step (z,k) + | k > h = return $ Done + | otherwise = return $ Yield (z:.Letter k) (z,k+1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamUp #-} + streamDown (ls:.Letter l) (hs:.Letter h) = flatten mk step Unknown $ streamDown ls hs + where mk z = return (z,h) + step (z,k) + | k < l = return $ Done + | otherwise = return $ Yield (z:.Letter k) (z,k-1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamDown #-} + +-- TODO temporary, because defaults dont inline + +instance IndexStream (Letter l) where + streamUp l h = map (\(Z:.k) -> k) $ streamUp (Z:.l) (Z:.h) + {-# Inline streamUp #-} + streamDown l h = map (\(Z:.k) -> k) $ streamDown (Z:.l) (Z:.h) + {-# Inline streamDown #-} + +{- +instance (Index sh, Show sh) => Shape (sh :. Letter z) where + rank (sh:._) = rank sh + 1 + zeroDim = zeroDim:.Letter 0 + unitDim = unitDim:.Letter 1 -- TODO does this one make sense? + intersectDim (sh1:.n1) (sh2:.n2) = intersectDim sh1 sh2 :. min n1 n2 + addDim (sh1:.Letter n1) (sh2:.Letter n2) = addDim sh1 sh2 :. Letter (n1+n2) -- TODO will not necessarily yield a valid Letter + size (sh1:.Letter n) = size sh1 * n + sizeIsValid (sh1:.Letter n) = sizeIsValid (sh1:.n) + toIndex (sh1:.Letter sh2) (sh1':.Letter sh2') = toIndex (sh1:.sh2) (sh1':.sh2') + fromIndex (ds:.Letter d) n = fromIndex ds (n `quotInt` d) :. Letter 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 :. Letter 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 :. Letter z) where + subDim (sh1:.Letter n1) (sh2:.Letter n2) = subDim sh1 sh2 :. Letter (n1-n2) + rangeList (sh1:.Letter n1) (sh2:.Letter n2) = [ sh:.Letter n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2)]] + rangeStream (fs:.Letter f) (ts:.Letter t) = VM.flatten mk step Unknown $ rangeStream fs ts where + mk sh = return (sh :. f) + step (sh :. k) + | k>t = return $ VM.Done + | otherwise = return $ VM.Yield (sh :. Letter k) (sh :. k +1) + {-# INLINE [1] mk #-} + {-# INLINE [1] step #-} + {-# INLINE rangeStream #-} +-} + 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,150 @@ + +{-# Language CPP #-} + +#if __GLASGOW_HASKELL__ < 710 +{-# LANGUAGE OverlappingInstances #-} +#endif + +-- | Convert between different nucleotide representations + +module Biobase.Primary.Nuc.Conversion where + +import qualified Data.Vector.Unboxed as VU + +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 #-} + + + +-- * 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,84 @@ + +module Biobase.Primary.Nuc.DNA where + +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.Primary.Bounds +import Biobase.Primary.Letter + + + +-- | DNA nucleotides. + +data DNA + +-- Single-character names for nucleotides. + +pattern A = Letter 0 :: Letter DNA +pattern C = Letter 1 :: Letter DNA +pattern G = Letter 2 :: Letter DNA +pattern T = Letter 3 :: Letter DNA +pattern N = Letter 4 :: Letter DNA + +instance Enum (Letter DNA) 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 + +acgt :: [Letter DNA] +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 #-} + +instance Show (Letter DNA) where + show c = [dnaChar c] + +instance Read (Letter DNA) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | otherwise = [(charDNA x, xs)] + +dnaSeq :: MkPrimary n DNA => n -> Primary DNA +dnaSeq = primary + +instance Bounded (Letter DNA) where + minBound = A + maxBound = N + +instance MkPrimary (VU.Vector Char) DNA where + primary = VU.map charDNA + +instance IsString [Letter DNA] 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,82 @@ + +module Biobase.Primary.Nuc.RNA where + +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.Primary.Bounds +import Biobase.Primary.Letter + + + +-- | RNA nucleotides. + +data RNA + +pattern A = Letter 0 :: Letter RNA +pattern C = Letter 1 :: Letter RNA +pattern G = Letter 2 :: Letter RNA +pattern U = Letter 3 :: Letter RNA +pattern N = Letter 4 :: Letter RNA + +instance Bounded (Letter RNA) where + minBound = A + maxBound = N + +instance Enum (Letter RNA) 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 + +acgu :: [Letter RNA] +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' +{-# INLINE rnaChar #-} + +instance Show (Letter RNA) where + show c = [rnaChar c] + +instance Read (Letter RNA) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | otherwise = [(charRNA x, xs)] + +rnaSeq :: MkPrimary n RNA => n -> Primary RNA +rnaSeq = primary + +instance MkPrimary (VU.Vector Char) RNA where + primary = VU.map charRNA + +instance IsString [Letter RNA] where + fromString = map charRNA + 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,82 @@ + +module Biobase.Primary.Nuc.XNA where + +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.Primary.Bounds +import Biobase.Primary.Letter + + + +-- | Combine both, RNA and DNA. + +data XNA + +pattern A = Letter 0 :: Letter XNA +pattern C = Letter 1 :: Letter XNA +pattern G = Letter 2 :: Letter XNA +pattern T = Letter 3 :: Letter XNA +pattern U = Letter 4 :: Letter XNA +pattern N = Letter 5 :: Letter XNA + +instance Bounded (Letter XNA) where + minBound = A + maxBound = N + +instance Enum (Letter XNA) 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 + +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) where + show c = [xnaChar c] + +instance Read (Letter XNA) where + readsPrec p [] = [] + readsPrec p (x:xs) + | x==' ' = readsPrec p xs + | otherwise = [(charXNA x, xs)] + +xnaSeq :: MkPrimary n XNA => n -> Primary XNA +xnaSeq = primary + +instance MkPrimary (VU.Vector Char) XNA where + primary = VU.map charXNA + +instance IsString [Letter XNA] where + fromString = map charXNA + 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,76 @@ + +-- | 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.Arrow ((***)) +import Data.ByteString.Char8 (ByteString,unpack) +import Data.FileEmbed (embedFile) +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.Primary.AA +import Biobase.Primary.Nuc +import Biobase.Primary.Letter + + + +-- | Using the codon table, create an amino acid sequence from a @DNA@ +-- sequence (encoded as 'Primary DNA'). Suffixed @seq@ as we deal with +-- sequences, not letters. + +dnaAAseq :: Primary DNA -> Primary AA +dnaAAseq = VU.fromList . go where + go (VU.length -> 0) = [] + go (VU.splitAt 3 -> (hs,ts)) = case M.lookup hs dnaAAmap of + Just aa -> aa : go ts + _ -> error $ "dnaAAseq: " ++ show (hs,ts) + +-- | Transform an amino acid sequence back into DNA. +-- +-- WARNING: This is lossy! + +aaDNAseq :: Primary AA -> Primary DNA +aaDNAseq = VU.concatMap go where + go aa = case M.lookup aa aaDNAmap of + Just codon -> codon + Nothing -> error $ "aaDNAseq" ++ show aa + + +-- * Embedded codon data + +-- | Lossy backtransformation. + +aaDNAmap :: M.Map (Letter AA) (Primary DNA) +aaDNAmap = M.fromList . map swap . M.assocs $ dnaAAmap +{-# NOINLINE aaDNAmap #-} + +dnaAAmap :: Map (Primary DNA) (Letter AA) +dnaAAmap = M.fromList . map (primary *** charAA) . M.assocs $ codonTable where +{-# NOINLINE dnaAAmap #-} + +codonTable :: Map String Char +codonTable = M.fromList . map (go . words) . lines . unpack $ codonListEmbedded where + go [cs,[c]] = (cs,c) + go e = error $ "codonTable:" ++ show e +{-# NOINLINE codonTable #-} + +-- | Raw codon table + +codonListEmbedded :: ByteString +codonListEmbedded = $(embedFile "sources/codontable") + 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 #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE TypeFamilies #-} --- | 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 Data.Array.Repa.Index -import 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 Data.Vector.Unboxed.Deriving -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' - -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' - -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 #-} - --- | simple cis/wc-wc basepairs - -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 #-} - -derivingUnbox "Edge" - [t| Edge -> Int |] [| unEdge |] [| Edge |] +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 -derivingUnbox "CTisomerism" - [t| CTisomerism -> Int |] [| unCT |] [| CT |] +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 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,312 @@ + +{-# LANGUAGE FunctionalDependencies #-} + +-- {-# LANGUAGE OverlappingInstances #-} + +-- | 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.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 +import Data.Tuple (swap) +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.Primary + + + +-- * 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 + +-- TODO Index instances! + +{- +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 #-} +-} + +-- | 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 + +-- TODO Index instances + +{- +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 #-} +-} + +-- | 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 = (Letter RNA,Letter RNA) + +-- | 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 + +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) + + + +-- * tuple-like selection +-- +-- the 'lens' library provides combinators that should make this +-- superfluous. + +-- | 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 #-} + +-- | simple cis/wc-wc basepairs + +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/Constraint.hs b/Biobase/Secondary/Constraint.hs --- a/Biobase/Secondary/Constraint.hs +++ b/Biobase/Secondary/Constraint.hs @@ -1,29 +1,23 @@ -{-# LANGUAGE PackageImports #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE StandaloneDeriving #-} +-- | Simple oldstyle RNAfold constraints. A constraint yields a bonus or +-- malus to energy. + module Biobase.Secondary.Constraint where -import Data.Array.Repa.Index -import Data.Array.Repa.Shape -import Data.Char (toLower) -import Data.Primitive.Types +import Data.Char (toLower) +import Data.Primitive.Types +import Prelude as P import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU -import Prelude as P -import Data.PrimitiveArray -import Data.PrimitiveArray.Zero +import Data.PrimitiveArray -import Biobase.Secondary.Diagrams +import Biobase.Secondary.Diagrams --- | We can create a constraint from different sources +-- | We can create a constraint from different sources. class MkConstraint a where mkConstraint :: a -> Constraint @@ -34,9 +28,11 @@ newtype Constraint = Constraint {unConstraint :: VU.Vector (Char,Int)} deriving (Show,Read,Eq) +bonusCC :: VU.Vector Char bonusCC = VU.fromList "()<>|" {-# NOINLINE bonusCC #-} +nobonusCC :: VU.Vector Char nobonusCC = VU.fromList ".x" {-# NOINLINE nobonusCC #-} @@ -52,9 +48,9 @@ -- TODO and again, we should parametrize over "Energy", "Score", etc (that is, -- Prim a) -bonusTable :: Double -> Double -> Constraint -> Unboxed DIM2 Double +bonusTable :: Double -> Double -> Constraint -> Unboxed (Z:.Int:.Int) Double bonusTable bonus malus (Constraint constraint) = arr where - arr = fromAssocs zeroDim (Z:.n:.n) 0 $ bonusBr ++ bonusAn ++ bonusBa ++ malusBr ++ malusAn ++ malusX + arr = fromAssocs (Z:.0:.0) (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 @@ -99,13 +95,6 @@ , 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 diff --git a/Biobase/Secondary/Diagrams.hs b/Biobase/Secondary/Diagrams.hs --- a/Biobase/Secondary/Diagrams.hs +++ b/Biobase/Secondary/Diagrams.hs @@ -1,48 +1,61 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE FlexibleInstances #-} -- | Types for RNA secondary structure. Types vary from the simplest array -- (D1Secondary) to rather complex ones. -{-# LANGUAGE RecordWildCards #-} - module Biobase.Secondary.Diagrams where - import Control.Applicative import Control.Arrow -import Data.List (sort,groupBy,sortBy) +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.Cereal +import GHC.Generics import qualified Data.Vector.Unboxed as VU import Text.Printf -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@. newtype D1Secondary = D1S {unD1S :: VU.Vector Int} - deriving (Read,Show,Eq) + deriving (Read,Show,Eq,Generic) +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 @@ -52,21 +65,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 @@ -75,8 +88,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 @@ -87,21 +100,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 @@ -124,7 +122,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 @@ -133,7 +131,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) @@ -149,16 +147,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 = unsafeDotBracket 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 = ')' @@ -173,22 +167,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 (`elem` "().") + +-- | Is constraint type structure, i.e. there can also be symbols present +-- that denote up- or downstream pairing. + +isConstraintStructure :: String -> Bool +isConstraintStructure = all (`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. -unsafeDotBracket :: [String] -> String -> [(Int,Int)] -unsafeDotBracket 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 [] = [] @@ -202,8 +242,8 @@ -- | Secondary structure parser with a notion of errors. We either return a -- @Right@ structure, including flags, or a @Left@ error. -dotBracket :: [String] -> String -> Either String ( [(Int,Int)] ) -dotBracket dict str = fmap (sort . concat) . sequence . map (f str) $ dict where +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 [] @@ -215,4 +255,13 @@ g a b c = fail $ printf "unspecified error: %s (dot-bracket: %s)" (show (a,b,c)) str f xs lr@(_:_:_:_) = fail $ printf "unsound dictionary: %s (dot-bracket: %s)" lr str f xs lr = fail $ 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 [] 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,27 @@ -{-# 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 (embedFile) +import Data.Function (on) +import Data.List 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 @@ -56,11 +55,11 @@ instance IsostericityLookup Pair 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 ((CWW==).baseT.fst) -- keep only cWW pairs (baseT-ype) . filter ((x `elem`).snd) -- select based on class $ M.assocs defaultIsostericityMap @@ -81,7 +80,7 @@ mkIsostericityList :: [[[String]]] -> [(ExtPair, [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 diff --git a/Biobase/Secondary/PseudoKnots.hs b/Biobase/Secondary/PseudoKnots.hs deleted file mode 100644 --- a/Biobase/Secondary/PseudoKnots.hs +++ /dev/null @@ -1,62 +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 IndexStream (z:.ViennaPair) where + streamUp (ls:.ViennaPair l) (hs:.ViennaPair h) = flatten mk step Unknown $ streamUp ls hs + where mk z = return (z,l) + step (z,k) + | k > h = return $ Done + | otherwise = return $ Yield (z:.ViennaPair k) (z,k+1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamUp #-} + streamDown (ls:.ViennaPair l) (hs:.ViennaPair h) = flatten mk step Unknown $ streamDown ls hs + where mk z = return (z,h) + step (z,k) + | k < l = return $ Done + | otherwise = return $ Yield (z:.ViennaPair k) (z,k-1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamDown #-} + +instance IndexStream ViennaPair where + streamUp l h = map (\(Z:.k) -> k) $ streamUp (Z:.l) (Z:.h) + {-# Inline streamUp #-} + streamDown l h = map (\(Z:.k) -> k) $ streamDown (Z:.l) (Z:.h) + {-# Inline streamDown #-} + + + +{- instance (Shape sh,Show sh) => Shape (sh :. ViennaPair) where rank (sh:._) = rank sh + 1 zeroDim = zeroDim:.ViennaPair 0 @@ -73,52 +112,67 @@ 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)]] + rangeStream (fs:.ViennaPair f) (ts:.ViennaPair t) = VM.flatten mk step Unknown $ rangeStream fs ts where + mk sh = return (sh :. f) + step (sh :. k) + | k>t = return $ VM.Done + | otherwise = return $ VM.Yield (sh :. ViennaPair k) (sh :. k +1) + {-# INLINE [1] mk #-} + {-# INLINE [1] step #-} + {-# INLINE rangeStream #-} +-} -(vpNP:vpCG:vpGC:vpGU:vpUG:vpAU:vpUA:vpNS:vpUndefined:_) = P.map ViennaPair [0..] +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 +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 #-} -isViennaPair :: Nuc -> Nuc -> Bool -isViennaPair a b = f a b where - f l r = l==nC && r==nG - || l==nG && r==nC - || l==nA && r==nU - || l==nU && r==nA - || l==nG && r==nU - || l==nU && r==nG - {-# INLINE f #-} +isViennaPair :: Letter RNA -> Letter RNA -> 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:.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) +viennaPairTable :: Unboxed (Z:.Letter RNA:.Letter RNA) ViennaPair +viennaPairTable = fromAssocs (Z:.N:.N) (Z:.U:.U) 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 #-} @@ -131,14 +185,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 @@ -159,23 +207,23 @@ -- | 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")] derivingUnbox "ViennaPair" [t| ViennaPair -> Int |] [| unViennaPair |] [| ViennaPair |] diff --git a/BiobaseXNA.cabal b/BiobaseXNA.cabal --- a/BiobaseXNA.cabal +++ b/BiobaseXNA.cabal @@ -1,16 +1,16 @@ name: BiobaseXNA -version: 0.8.3.0 +version: 0.9.1.0 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-2014 +maintainer: choener@bioinf.uni-leipzig.de +homepage: http://www.bioinf.uni-leipzig.de/~choener/ +copyright: Christian Hoener zu Siederdissen, 2011 - 2015 category: Bioinformatics synopsis: Efficient RNA/DNA representations license: GPL-3 license-file: LICENSE build-type: Simple stability: experimental -cabal-version: >= 1.6.0 +cabal-version: >= 1.10.0 description: This is a base library for bioinformatics with emphasis on RNA and DNA primary structure as well as amino acid sequences. @@ -35,45 +35,88 @@ sources/isostericity-matrices.csv sources/isostericity-detailed.csv sources/iupac-nucleotides - changelog + sources/codontable + changelog.md library - build-depends: - base >3 && <5, - bytestring >= 0.10 , - containers >= 0.4 , - csv >= 0.1.2 , - file-embed >= 0.0.4.7 , - mtl >= 2.1 , - primitive >= 0.5 , - PrimitiveArray >= 0.5.4 , - repa >= 3.2 , - text >= 0.11 , - tuple >= 0.2 , - vector >= 0.10 , - vector-th-unbox >= 0.2 + build-depends: base >= 4.7 && < 4.9 + , aeson == 0.8.* + , bimaps == 0.0.0.* + , binary == 0.7.* + , bytes == 0.15.* + , bytestring == 0.10.* + , cereal == 0.4.* + , cereal-vector == 0.2.* + , containers == 0.5.* + , csv == 0.1.* + , file-embed == 0.0.8.* + , hashable == 1.2.* + , lens == 4.* + , primitive >= 0.5 && < 0.7 + , PrimitiveArray == 0.6.* + , split == 0.2.* + , text == 1.* + , tuple == 0.3.* + , vector == 0.10.* + , vector-binary-instances == 0.2.* + , vector-th-unbox == 0.2.* exposed-modules: - Biobase.AAseq Biobase.Primary - Biobase.Primary.Bounds + Biobase.Primary.AA Biobase.Primary.Hashed Biobase.Primary.IUPAC + Biobase.Primary.Letter + Biobase.Primary.Nuc + Biobase.Primary.Nuc.Conversion + Biobase.Primary.Nuc.DNA + Biobase.Primary.Nuc.RNA + Biobase.Primary.Nuc.XNA + Biobase.Primary.Trans Biobase.Secondary + Biobase.Secondary.Basepair Biobase.Secondary.Constraint Biobase.Secondary.Diagrams Biobase.Secondary.Isostericity - Biobase.Secondary.PseudoKnots + Biobase.Secondary.Pseudoknots + Biobase.Secondary.Structure Biobase.Secondary.Vienna + default-extensions: BangPatterns + , DeriveGeneric + , EmptyDataDecls + , FlexibleContexts + , FlexibleInstances + , GeneralizedNewtypeDeriving + , LambdaCase + , MultiParamTypeClasses + , PatternSynonyms + , ScopedTypeVariables + , TemplateHaskell + , TypeFamilies + , TypeOperators + , UndecidableInstances + , ViewPatterns + + default-language: Haskell2010 + ghc-options: -O2 -funbox-strict-fields executable SubOptDistance - build-depends: - cmdargs >= 0.10 + build-depends: base + , cmdargs == 0.10.* + , BiobaseXNA main-is: SubOptDistance.hs + hs-source-dirs: + src + default-language: + Haskell2010 + default-extensions: DeriveDataTypeable + , NoMonomorphismRestriction + , RecordWildCards + , ScopedTypeVariables ghc-options: -O2 diff --git a/SubOptDistance.hs b/SubOptDistance.hs deleted file mode 100644 --- a/SubOptDistance.hs +++ /dev/null @@ -1,53 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE DeriveDataTypeable #-} -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE NoMonomorphismRestriction #-} - -module Main where - -import Control.Arrow -import Data.Char (isSpace) -import Data.Either (either) -import Data.List -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 a bunch of RNAsubopt lines 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 - let xs = map (\(s,e) -> (s,either error id $ dotBracket ["()"] s,read e)) $ map (break isSpace) xs' - let q = either error id $ dotBracket ["()"] structure - let (d,x,e) = getMinimalDistance aq qa q xs - let distance :: Double = e - ((/100) . read $ words sqn !! 1) - printf "%4d %s %8.2f %8.2f\n" d x e distance - -getMinimalDistance :: Bool -> Bool -> [PairIdx] -> [(String,[PairIdx],Double)] -> (Int,String,Double) -getMinimalDistance aq qa q xs = g $ minimumBy (comparing f) xs where - g (a,b,c) = ( fst $ f (a,b,c) , a, c) - f (_,a,e) = (length $ aqs ++ qas , e) where - aqs = if aq then a \\ q else [] - qas = if qa then q \\ a else [] diff --git a/changelog b/changelog deleted file mode 100644 --- a/changelog +++ /dev/null @@ -1,37 +0,0 @@ -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/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,58 @@ +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/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 +