BiobaseXNA 0.8.3.0 → 0.9.1.0
raw patch · 28 files changed
+1885/−1244 lines, 28 filesdep +BiobaseXNAdep +aesondep +bimapsdep −mtldep −repadep ~PrimitiveArraydep ~basedep ~bytestring
Dependencies added: BiobaseXNA, aeson, bimaps, binary, bytes, cereal, cereal-vector, hashable, lens, split, vector-binary-instances
Dependencies removed: mtl, repa
Dependency ranges changed: PrimitiveArray, base, bytestring, cmdargs, containers, csv, file-embed, primitive, text, tuple, vector, vector-th-unbox
Files
- Biobase/AAseq.hs +0/−254
- Biobase/Primary.hs +16/−200
- Biobase/Primary/AA.hs +152/−0
- Biobase/Primary/Bounds.hs +0/−52
- Biobase/Primary/Hashed.hs +23/−29
- Biobase/Primary/IUPAC.hs +127/−24
- Biobase/Primary/Letter.hs +158/−0
- Biobase/Primary/Nuc.hs +28/−0
- Biobase/Primary/Nuc/Conversion.hs +150/−0
- Biobase/Primary/Nuc/DNA.hs +84/−0
- Biobase/Primary/Nuc/RNA.hs +82/−0
- Biobase/Primary/Nuc/XNA.hs +82/−0
- Biobase/Primary/Trans.hs +76/−0
- Biobase/Secondary.hs +16/−345
- Biobase/Secondary/Basepair.hs +312/−0
- Biobase/Secondary/Constraint.hs +13/−24
- Biobase/Secondary/Diagrams.hs +101/−52
- Biobase/Secondary/Isostericity.hs +13/−14
- Biobase/Secondary/PseudoKnots.hs +0/−62
- Biobase/Secondary/Pseudoknots.hs +60/−0
- Biobase/Secondary/Structure.hs +36/−0
- Biobase/Secondary/Vienna.hs +121/−73
- BiobaseXNA.cabal +68/−25
- SubOptDistance.hs +0/−53
- changelog +0/−37
- changelog.md +58/−0
- sources/codontable +64/−0
- src/SubOptDistance.hs +45/−0
− Biobase/AAseq.hs
@@ -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<n2 && inShapeRange sh1 sh2 idx- listOfShape (sh:.AA n) = n : listOfShape sh- shapeOfList xx = case xx of- [] -> 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 |]-
Biobase/Primary.hs view
@@ -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<n2 && inShapeRange sh1 sh2 idx- listOfShape (sh:.Nuc n) = n : listOfShape sh- shapeOfList xx = case xx of- [] -> 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
+ Biobase/Primary/AA.hs view
@@ -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+
− Biobase/Primary/Bounds.hs
@@ -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)
Biobase/Primary/Hashed.hs view
@@ -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 #-}
Biobase/Primary/IUPAC.hs view
@@ -1,47 +1,150 @@-{-# LANGUAGE TemplateHaskell #-} -- | Degenerate base symbol representation. We use the same conventions as in -- <<https://en.wikipedia.org/wiki/Nucleic_acid_notation>> 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.
+ Biobase/Primary/Letter.hs view
@@ -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<n2 && inShapeRange sh1 sh2 idx+ listOfShape (sh:.Letter n) = n : listOfShape sh+ shapeOfList xx = case xx of+ [] -> 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 #-}+-}+
+ Biobase/Primary/Nuc.hs view
@@ -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+
+ Biobase/Primary/Nuc/Conversion.hs view
@@ -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+
+ Biobase/Primary/Nuc/DNA.hs view
@@ -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+
+ Biobase/Primary/Nuc/RNA.hs view
@@ -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+
+ Biobase/Primary/Nuc/XNA.hs view
@@ -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+
+ Biobase/Primary/Trans.hs view
@@ -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")+
Biobase/Secondary.hs view
@@ -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<n2 && inShapeRange sh1 sh2 idx- listOfShape (sh:.Edge n) = n : listOfShape sh- shapeOfList xx = case xx of- [] -> 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<n2 && inShapeRange sh1 sh2 idx- listOfShape (sh:.CT n) = n : listOfShape sh- shapeOfList xx = case xx of- [] -> 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
+ Biobase/Secondary/Basepair.hs view
@@ -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<n2 && inShapeRange sh1 sh2 idx+ listOfShape (sh:.Edge n) = n : listOfShape sh+ shapeOfList xx = case xx of+ [] -> 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<n2 && inShapeRange sh1 sh2 idx+ listOfShape (sh:.CT n) = n : listOfShape sh+ shapeOfList xx = case xx of+ [] -> 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 #-}+
Biobase/Secondary/Constraint.hs view
@@ -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
Biobase/Secondary/Diagrams.hs view
@@ -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 []
Biobase/Secondary/Isostericity.hs view
@@ -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
− Biobase/Secondary/PseudoKnots.hs
@@ -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<k&&k<j&&j<l || k<i&&i<l&&l<j) ps--instance RemovePseudoKnots [PairIdx] where- removeByCounting = VU.toList . removeByCounting . VU.fromList---- | Remove pseudoknotted pairs from extended RNA secondary structures.--instance RemovePseudoKnots (VU.Vector ExtPairIdx) 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<k&&k<j&&j<l || k<i&&i<l&&l<j) ps--instance RemovePseudoKnots [ExtPairIdx] where- removeByCounting = VU.toList . removeByCounting . VU.fromList-
+ Biobase/Secondary/Pseudoknots.hs view
@@ -0,0 +1,60 @@++-- | 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 Data.List+import qualified Data.Vector.Unboxed as VU++import Biobase.Secondary.Basepair++++-- | 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<k&&k<j&&j<l || k<i&&i<l&&l<j) ps++instance RemovePseudoKnots [PairIdx] where+ removeByCounting = VU.toList . removeByCounting . VU.fromList++-- | Remove pseudoknotted pairs from extended RNA secondary structures.++instance RemovePseudoKnots (VU.Vector ExtPairIdx) 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<k&&k<j&&j<l || k<i&&i<l&&l<j) ps++instance RemovePseudoKnots [ExtPairIdx] where+ removeByCounting = VU.toList . removeByCounting . VU.fromList+
+ Biobase/Secondary/Structure.hs view
@@ -0,0 +1,36 @@++-- | A secondary structure, with sequence, Vienna compatible canonical+-- secondary structure, extended structure, and additional information.+--+-- This is the structure that will be returned by prediction algorithms in+-- the future.+--+-- TODO we will need ex- and import functions to a number of standard+-- formats. There is an open feature request to export to something that+-- resembles FASTA with additional information.++module Biobase.Secondary.Structure where++import Data.Map.Strict (Map)+import Data.Text (Text)+import qualified Data.Text as T++import Biobase.Secondary.Diagrams++++-- | A sequence, complete with secondary structure. While this structure is+-- rather RNA-centric, there is nothing that prohibits us from using this+-- for DNA.+--+-- TODO Generics, Cereal, Binary, Aeson instances++data SecondaryStructure = SS+ { _ssSeq :: !Text -- ^ sequence; we use 'Text' whenever possible+ , _ssVienna :: !D1Secondary -- ^ canonical Vienna secondary structure+ , _ssExt :: !D2Secondary -- ^ extended secondary structure+ , _ssViennaE :: Maybe () -- ^ TODO will be the energy, measured or predicted+ , _ssAux :: Map Text Text -- ^ any auxiliary info in key/value format+ }+ deriving (Eq,Show,Read)+
Biobase/Secondary/Vienna.hs view
@@ -1,43 +1,82 @@-{-# LANGUAGE PackageImports #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE PatternGuards #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE StandaloneDeriving #-}-{-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE TypeFamilies #-} -- | Encoding of Watson-Crick and Wobble Pairs in the Vienna RNA package style. module Biobase.Secondary.Vienna where -import Data.Array.Repa.Index-import Data.Array.Repa.Shape+import Data.Aeson+import Data.Binary import Data.Ix import Data.Primitive.Types+import Data.Serialize (Serialize(..)) import Data.Tuple (swap)+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 Prelude as P+import GHC.Generics (Generic)+import Prelude hiding (map) import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Unboxed as VU+import qualified Prelude as P -import Data.Array.Repa.ExtShape-import Data.PrimitiveArray as PA-import Data.PrimitiveArray.Zero as PA+import Data.PrimitiveArray hiding (Complement(..),map) -import Biobase.Primary-import Biobase.Primary.Bounds+import Biobase.Primary.Letter+import Biobase.Primary.Nuc+import Biobase.Primary.Nuc.RNA -- | Use machine Ints internally newtype ViennaPair = ViennaPair { unViennaPair :: Int }- deriving (Eq,Ord,Ix)+ deriving (Eq,Ord,Generic,Ix) +instance Binary (ViennaPair)+instance Serialize (ViennaPair)+instance FromJSON (ViennaPair)+instance ToJSON (ViennaPair)++instance Index ViennaPair where+ linearIndex _ _ (ViennaPair p) = p+ {-# Inline linearIndex #-}+ smallestLinearIndex _ = error "still needed?"+ {-# Inline smallestLinearIndex #-}+ largestLinearIndex (ViennaPair p) = p+ {-# Inline largestLinearIndex #-}+ size _ (ViennaPair h) = h+1+ {-# Inline size #-}+ inBounds (ViennaPair l) (ViennaPair h) (ViennaPair p) = l <= p && p <= h+ {-# Inline inBounds #-}++instance IndexStream z => 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 |]
BiobaseXNA.cabal view
@@ -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
− SubOptDistance.hs
@@ -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 []
− changelog
@@ -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
+ changelog.md view
@@ -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
+ sources/codontable view
@@ -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
+ src/SubOptDistance.hs view
@@ -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+