diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License + +Copyright (c) 2008-2011 Nicholas Ingolia + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/seqloc.cabal b/seqloc.cabal new file mode 100644 --- /dev/null +++ b/seqloc.cabal @@ -0,0 +1,43 @@ +Name: seqloc +Version: 0.1 +Cabal-Version: >= 1.2 +Synopsis: Handle sequence locations for bioinformatics +Description: Handle sequence locations for bioinformatics +Homepage: http://www.ingolia-lab.org/software/seqloc/ +License: BSD3 +License-file: LICENSE +Author: Nick Ingolia +Maintainer: nick@ingolia.org +Build-Type: Simple +Cabal-Version: >= 1.2 + +Category: Bioinformatics + +Flag Tests + Description: Build test program + Default: False + +Library + Exposed-modules: Bio.SeqLoc.Strand, Bio.SeqLoc.Position, + Bio.SeqLoc.SeqLike, Bio.SeqLoc.Location, + Bio.SeqLoc.SpliceLocation, + Bio.SeqLoc.OnSeq, + Bio.SeqLoc.LocRepr, + Bio.SeqLoc.Transcript + Build-depends: base >= 4.3 && < 5, bytestring, haskell98, attoparsec >= 0.8.5 + Hs-Source-Dirs: src + Ghc-options: -Wall + +Executable test-seqloc + if !flag(Tests) + Buildable: False + Main-is: TestMain.hs + Other-modules: Bio.SeqLoc.Strand, Bio.SeqLoc.Position, + Bio.SeqLoc.SeqLike, Bio.SeqLoc.Location, + Bio.SeqLoc.SpliceLocation, + Bio.SeqLoc.OnSeq, + Bio.SeqLoc.LocRepr + Build-depends: base >= 4.3 && < 5, bytestring, haskell98, attoparsec >= 0.8.5, + QuickCheck, random + Hs-Source-Dirs: src, test + Ghc-options: -Wall diff --git a/src/Bio/SeqLoc/LocRepr.hs b/src/Bio/SeqLoc/LocRepr.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/LocRepr.hs @@ -0,0 +1,25 @@ +module Bio.SeqLoc.LocRepr + ( LocRepr(..), reprStr + , unreprMaybe, unreprEither, unreprErr + ) + where + +import qualified Data.ByteString.Char8 as BS + +import qualified Data.Attoparsec.Zepto as ZP + +class LocRepr l where + repr :: l -> BS.ByteString + unrepr :: ZP.Parser l + +reprStr :: (LocRepr l) => l -> String +reprStr = BS.unpack . repr + +unreprMaybe :: (LocRepr l) => BS.ByteString -> Maybe l +unreprMaybe = either (const Nothing) Just . ZP.parse unrepr + +unreprEither :: (LocRepr l) => BS.ByteString -> Either String l +unreprEither = ZP.parse unrepr + +unreprErr :: (LocRepr l) => BS.ByteString -> l +unreprErr = either error id . ZP.parse unrepr diff --git a/src/Bio/SeqLoc/Location.hs b/src/Bio/SeqLoc/Location.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/Location.hs @@ -0,0 +1,237 @@ +{-# LANGUAGE TypeFamilies, FlexibleContexts #-} +{-| Data type for a sequence location consiting of a contiguous range +of positions on the sequence. + +Throughout, /sequence position/ refers to a 'Pos.Pos' which includes a +strand, as opposed to an /offset/, which refers to a 'Pos.Offset' with +no strand. + + -} + +module Bio.SeqLoc.Location + ( -- * Sequence locations + Location(..), overlaps + -- * Contiguous sequence locations + , ContigLoc, offset5, fromStartEnd, fromPosLen, fromBoundsStrand + -- * Transforming locations + , slide + ) + where + +import Prelude hiding (length) + +import Control.Applicative +import Control.Monad +import qualified Data.ByteString.Char8 as BS + +import qualified Data.Attoparsec.Zepto as ZP + +import Bio.SeqLoc.LocRepr +import qualified Bio.SeqLoc.Position as Pos +import Bio.SeqLoc.Strand +import qualified Bio.SeqLoc.SeqLike as SeqLike + +class Location l where + strand :: l -> Strand + length :: l -> Pos.Offset + + -- | The bounds of a sequence location. This is a pair consisting + -- of the lowest and highest sequence offsets covered by the region. + -- The bounds ignore the strand of the sequence location, and the + -- first element of the pair will always be lower than the second. + bounds :: l -> (Pos.Offset, Pos.Offset) + + -- | Sequence position of the start of the location. This is the 5' + -- end on the location strand, which will have a higher offset than + -- 'endPos' if the location is on the 'RevCompl' strand. + startPos :: l -> Pos.Pos + + -- | Sequence position of the end of the location, as described in + -- 'startPos'. + endPos :: l -> Pos.Pos + + -- | Extract 'Just' the nucleotide 'SeqLike' for the sequence + -- location, or 'Nothing' if f any part of the location lies outside + -- the bounds of the sequence. + seqData :: (SeqLike.SeqLike s, Stranded s) => s -> l -> Maybe s + + -- | As 'seqData', extract the nucleotide subsequence for the + -- location, but any positions in the location lying outside the + -- bounds of the sequence are returned as @N@. + seqDataPad :: (SeqLike.SeqLike s, Stranded s) => s -> l -> s + + -- | Given a sequence position and a sequence location relative to + -- the same sequence, compute a new position representing the + -- original position relative to the subsequence defined by the + -- location. If the sequence position lies outside of the sequence + -- location, @Nothing@ is returned; thus, the offset of the new + -- position will always be in the range @[0, length l - 1]@. + posInto :: Pos.Pos -> l -> Maybe Pos.Pos + + -- | Given a sequence location and a sequence position within that + -- location, compute a new position representing the original + -- position relative to the outer sequence. If the sequence + -- position lies outside the location, @Nothing@ is returned. + -- + -- This function inverts 'posInto' when the sequence position lies + -- within the position is actually within the location. + posOutof :: Pos.Pos -> l -> Maybe Pos.Pos + + -- | For an enclosing location and a sublocation in the same + -- coordinate system, find the image of the sublocation relative to + -- the enclosing location. For example, if the enclosing location + -- spans (100, 150) and the sublocation is (110, 120) then + -- 'clocInto' will be (10, 20). + clocInto :: ContigLoc -> l -> Maybe ContigLoc + + -- | Returns a sequence location produced by finding the inverse + -- image of a sublocation, with coordinates given relative to an + -- enclosing location, in the coordinate system of the enclosing + -- location. For example, if the enclosing location spans (100, + -- 150) and the sublocation is (10, 20) then 'clocOutof' will be + -- (110, 120). + clocOutof :: ContigLoc -> l -> Maybe l + + -- | Returns a sequence location produced by extending the original + -- location on each end, based on a pair of (/5\' extension/, /3\' + -- extension/). The 5\' extension is applied to the 5\' end of the + -- location on the location strand; if the location is on the + -- 'RevCompl' strand, the 5\' end will have a higher offset than the + -- 3\' end and this offset will increase by the amount of the 5\' + -- extension. Similarly, the 3\' extension is applied to the 3\' + -- end of the location. + extend :: (Pos.Offset, Pos.Offset) -> l -> l + + -- | Returns @True@ when a sequence offset lies within a sequence + -- location on the same sequence + offsetWithin :: Pos.Offset -> l -> Bool + + -- | Returns @True@ when a sequence position lies within a sequence + -- location on the same sequence, and occupies the same strand. + posWithin :: Pos.Pos -> l -> Bool + + -- | Returns @True@ when two sequence locations overlap at any + -- position. + contigOverlaps :: ContigLoc -> l -> Bool + + -- | Contigs that comprise the location + toContigs :: l -> [ContigLoc] + + +overlaps :: (Location l1, Location l2) => l1 -> l2 -> Bool +overlaps l1 = any (\c2 -> contigOverlaps c2 l1) . toContigs + +-- | Contiguous sequence location defined by a span of sequence +-- positions, lying on a specific strand of the sequence. +data ContigLoc = ContigLoc { offset5 :: !Pos.Offset -- ^ The offset of the 5\' end of the location, as a 0-based index + , clocLength :: !Pos.Offset -- ^ The length of the location + , clocStrand :: !Strand -- ^ The strand of the location + } deriving (Eq, Ord, Show) + +instance Stranded ContigLoc where + revCompl (ContigLoc seq5 len str) = ContigLoc seq5 len $ revCompl str + +to :: BS.ByteString +to = BS.pack "to" + +instance LocRepr ContigLoc where + repr cloc = let (seq5, seq3) = bounds cloc + in BS.concat [ repr seq5, to, repr seq3, repr . strand $ cloc ] + unrepr = fromBoundsStrand <$> unrepr <* ZP.string to <*> unrepr <*> unrepr + +instance Location ContigLoc where + strand = clocStrand + length = clocLength + seqData sequ (ContigLoc seq5 len str) = liftM (stranded str) . (SeqLike.subseq seq5 len) $ sequ + seqDataPad sequ (ContigLoc seq5 len str) = (stranded str) . (SeqLike.subseqPad seq5 len) $ sequ + posInto = clocPosInto + posOutof = clocPosOutof + bounds (ContigLoc seq5 len _) = (seq5, seq5 + len - 1) + startPos (ContigLoc seq5 len str) + = case str of + Fwd -> Pos.Pos seq5 str + RevCompl -> Pos.Pos (seq5 + len - 1) str + endPos (ContigLoc seq5 len str) + = case str of + Fwd -> Pos.Pos (seq5 + len - 1) str + RevCompl -> Pos.Pos seq5 str + clocInto = clocClocInto + clocOutof = clocClocOutof + extend = clocExtend + offsetWithin off (ContigLoc seq5 len _) + = (off >= seq5) && (off < seq5 + len) + posWithin (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand) + = (pos >= seq5) && (pos < seq5 + len) && (cStrand == pStrand) + contigOverlaps = clocOverlaps + toContigs = (: []) + +-- | Create a sequence location between 0-based starting and ending +-- bounds with a specified strand. +fromBoundsStrand :: Pos.Offset -> Pos.Offset -> Strand -> ContigLoc +fromBoundsStrand seq5 seq3 _ | seq3 < seq5 = error "Bio.SeqLoc.Location.fromBoundsStrand: seq3 < seq5" +fromBoundsStrand seq5 seq3 str = ContigLoc seq5 (1 + seq3 - seq5) str + +-- | Create a sequence location lying between 0-based starting and +-- ending offsets. When @start < end@, the location +-- be on the forward strand, otherwise it will be on the +-- reverse complement strand. +fromStartEnd :: Pos.Offset -> Pos.Offset -> ContigLoc +fromStartEnd start end + | start < end = ContigLoc start (1 + end - start) Fwd + | otherwise = ContigLoc end (1 + start - end) RevCompl + +-- | Create a sequence location from the sequence position of the +-- start of the location and the length of the position. The strand +-- of the location, and the direction it extends from the starting +-- position, are determined by the strand of the starting position. +fromPosLen :: Pos.Pos -> Pos.Offset -> ContigLoc +fromPosLen _ len | len < 0 = error "Bio.SeqLoc.Location.fromPosLen: len < 0" +fromPosLen (Pos.Pos off5 Fwd) len = ContigLoc off5 len Fwd +fromPosLen (Pos.Pos off3 RevCompl) len = ContigLoc (off3 - (len - 1)) len RevCompl + +-- | Returns a location resulting from sliding the original location +-- along the sequence by a specified offset. A positive offset will +-- move the location away from the 5\' end of the forward stand of the +-- sequence regardless of the strand of the location itself. Thus, +-- +-- > slide (revCompl cloc) off == revCompl (slide cloc off) +slide :: Pos.Offset -> ContigLoc -> ContigLoc +slide dpos (ContigLoc seq5 len str) = ContigLoc (seq5 + dpos) len str + +clocPosInto :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos +clocPosInto (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand) + | pos < seq5 || pos >= seq5 + len = Nothing + | otherwise = Just $ case cStrand of + Fwd -> Pos.Pos (pos - seq5) pStrand + RevCompl -> Pos.Pos (seq5 + len - (pos + 1)) (revCompl pStrand) + +clocPosOutof :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos +clocPosOutof (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand) + | pos < 0 || pos >= len = Nothing + | otherwise = Just $ case cStrand of + Fwd -> Pos.Pos (pos + seq5) pStrand + RevCompl -> Pos.Pos (seq5 + len - (pos + 1)) (revCompl pStrand) + +clocExtend :: (Pos.Offset, Pos.Offset) -> ContigLoc -> ContigLoc +clocExtend (ext5, ext3) (ContigLoc seq5 len str) + = case str of + Fwd -> ContigLoc (seq5 - ext5) (len + ext5 + ext3) str + RevCompl -> ContigLoc (seq5 - ext3) (len + ext5 + ext3) str + +clocClocInto :: ContigLoc -> ContigLoc -> Maybe ContigLoc +clocClocInto subcloc tocloc + = case (posInto (startPos subcloc) tocloc, posInto (endPos subcloc) tocloc) of + (Just start, Just _) -> Just (fromPosLen start (length subcloc)) + _ -> Nothing + +clocClocOutof :: ContigLoc -> ContigLoc -> Maybe ContigLoc +clocClocOutof subcloc fromcloc + = case (posOutof (startPos subcloc) fromcloc, posOutof (endPos subcloc) fromcloc) of + (Just start, Just _) -> Just (fromPosLen start (length subcloc)) + _ -> Nothing + +clocOverlaps :: ContigLoc -> ContigLoc -> Bool +clocOverlaps contig1 contig2 + = case (bounds contig1, bounds contig2) of + ((low1, high1),(low2, high2)) -> (strand contig1 == strand contig2) + && (low1 <= high2) && (low2 <= high1) diff --git a/src/Bio/SeqLoc/OnSeq.hs b/src/Bio/SeqLoc/OnSeq.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/OnSeq.hs @@ -0,0 +1,71 @@ +{-# LANGUAGE FlexibleInstances, TypeFamilies, FlexibleContexts, GeneralizedNewtypeDeriving #-} + +{-| + +Data types for sequence locations and sequence positions associated +with specific, named sequences. + +-} + +module Bio.SeqLoc.OnSeq ( + SeqName(..) + + , OnSeq(..) + + -- * Positions on named sequences + , SeqOffset, SeqPos + + -- * Contiguous location spans on named sequences + , ContigSeqLoc + + -- * Arbitrary location spans on named sequences + , SpliceSeqLoc + + , andSameSeq + + ) + where + +import Control.Applicative +import qualified Data.ByteString.Char8 as BS +import Data.ByteString.Internal (c2w) +import Data.String + +import qualified Data.Attoparsec.Zepto as ZP + +import Bio.SeqLoc.LocRepr +import qualified Bio.SeqLoc.Location as Loc +import qualified Bio.SeqLoc.Position as Pos +import qualified Bio.SeqLoc.SpliceLocation as SpLoc +import Bio.SeqLoc.Strand + +newtype SeqName = SeqName { unSeqName :: BS.ByteString } deriving (Show, Read, Eq, Ord, IsString) + +data OnSeq s = OnSeq { onSeqName :: !SeqName, unOnSeq :: !s } deriving (Show, Read, Eq, Ord) + +at :: BS.ByteString +at = BS.singleton '@' + +instance Stranded s => Stranded (OnSeq s) where + revCompl (OnSeq name obj) = OnSeq name (revCompl obj) + +instance LocRepr s => LocRepr (OnSeq s) where + repr (OnSeq name obj) = BS.concat [ unSeqName name, at, repr obj ] + unrepr = OnSeq <$> (SeqName <$> ZP.takeWhile (/= c2w '@')) <* ZP.string at <*> unrepr + +type SeqOffset = OnSeq Pos.Offset + +-- | A position on a named sequence +type SeqPos = OnSeq Pos.Pos + +-- | A location consisting of a contiguous span of positions on a +-- named sequence. +type ContigSeqLoc = OnSeq Loc.ContigLoc + +-- | A general location, consisting of spans of sequence positions on +-- a specific, named sequence. +type SpliceSeqLoc = OnSeq SpLoc.SpliceLoc + +andSameSeq :: (a -> b -> Bool) -> OnSeq a -> OnSeq b -> Bool +andSameSeq p (OnSeq n1 x) (OnSeq n2 y) | n1 == n2 = p x y + | otherwise = False diff --git a/src/Bio/SeqLoc/Position.hs b/src/Bio/SeqLoc/Position.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/Position.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE FlexibleContexts, GeneralizedNewtypeDeriving, OverloadedStrings #-} + +{-| Data type for a sequence position. + +Zero-based 'Offset'indices are used throughout, to facilitate direct +use of indexing functions on 'SeqData'. + +-} + +module Bio.SeqLoc.Position ( + -- * Sequence positions + Offset(..) + , Pos(..) + + -- * Manipulating positions + , slide + + -- * Extracting sequences + , atPos + ) + where + +import Control.Applicative +import Control.Monad (liftM) +import qualified Data.ByteString as BSW +import qualified Data.ByteString.Char8 as BS + +import qualified Data.Attoparsec.Char8 as AP (isDigit_w8) +import qualified Data.Attoparsec.Zepto as ZP + +import Bio.SeqLoc.LocRepr +import Bio.SeqLoc.SeqLike +import Bio.SeqLoc.Strand + +-- | Unstranded offset in a sequence +newtype Offset = Offset { unOffset :: Int } deriving (Eq, Ord, Show, Read, Num, Real, Enum, Integral) + +instance LocRepr Offset where + repr = BS.pack . show . unOffset + unrepr = (negate <$> (ZP.string "-" *> decimal)) <|> (ZP.string "+" *> decimal) <|> decimal + where decimal = Offset . BSW.foldl' step 0 <$> ZP.takeWhile AP.isDigit_w8 + step a w = a * 10 + fromIntegral (w - 48) + +-- | Stranded position in a sequence +data Pos = Pos { offset :: !Offset -- ^ 0-based index of the position + , strand :: !Strand -- ^ Strand of the position + } + deriving (Eq, Ord, Show, Read) + +instance Stranded Pos where + revCompl (Pos off str) = Pos off (revCompl str) + +instance LocRepr Pos where + repr (Pos off str) = BS.concat [ repr off, repr str ] + unrepr = Pos <$> unrepr <*> unrepr + +-- | Returns a position resulting from sliding the original position +-- along the sequence by a specified offset. A positive offset will +-- move the position away from the 5\' end of the forward stand of the +-- sequence regardless of the strand of the position itself. Thus, +-- +-- > slide (revCompl pos) off == revCompl (slide pos off) +slide :: Pos -> Offset -> Pos +slide (Pos off str) doff = Pos (off + doff) str + +-- | Extract 'Just' the item at a specific sequence position, or +-- 'Nothing' if the position lies outside the bounds of the sequence. +atPos :: (SeqLike s) => s -> Pos -> Maybe Char +atPos sequ (Pos off str) = liftM (stranded str) . ntAt sequ $ off + +{-# SPECIALIZE atPos :: String -> Pos -> Maybe Char #-} +{-# SPECIALIZE atPos :: BS.ByteString -> Pos -> Maybe Char #-} diff --git a/src/Bio/SeqLoc/SeqLike.hs b/src/Bio/SeqLoc/SeqLike.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/SeqLike.hs @@ -0,0 +1,127 @@ +{-# LANGUAGE FlexibleInstances #-} + +module Bio.SeqLoc.SeqLike + ( SeqLike(..) + ) + where + +import Prelude hiding (length, concat) + +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as LBS +import qualified Data.List as L +import Data.Maybe + +class SeqLike s where + -- | Length of sequence data + length :: (Integral n) => s -> n + + -- | 'Just' the nucleotide at a specified sequence data offset, + -- given in 0-based coordinates, or 'Nothing' if the offset is + -- beyond the bounds of the data + ntAt :: (Integral n) => s -> n -> Maybe Char + + -- | 'Just' the nucleotides in subsequence of the sequence data, or + -- 'Nothing' if the region extends beyond the bounds of the + -- sequence. + subseq :: (Integral n, Integral m) + => n -- ^ Starting position in 0-based coordinates + -> m -- ^ Length + -> s -- ^ Sequence data + -> Maybe s + + -- | Nucleotides in a subsequence of the sequence data, padded with + -- @N@ when the region extends beyond the bounds of the sequence. + subseqPad :: (Integral n, Integral m) + => n -- ^ Starting position in 0-based coordinates + -> m -- ^ Length + -> s -- ^ Sequence data + -> s + + concat :: [s] -> s + +instance SeqLike [Char] where + length = L.genericLength + ntAt l p | p < 0 = Nothing + | otherwise = listToMaybe . drop (fromIntegral p) $ l + subseq = listSubseq + subseqPad = listSubseqPad + concat = L.concat + +listSubseq :: (Integral n, Integral m) => n -> m -> [Char] -> Maybe [Char] +listSubseq start len l | start < 0 = Nothing + | otherwise = case take (fromIntegral len) . drop (fromIntegral start) $ l of + ss | L.genericLength ss == len -> Just ss + | otherwise -> Nothing + +listSubseqPad :: (Integral n, Integral m) => n -> m -> [Char] -> [Char] +listSubseqPad istart ilen sequ + | start + len <= 0 = replicate len 'N' + | start < 0 = replicate (negate start) 'N' ++ takePadded (len + start) sequ + | otherwise = takePadded len . drop start $ sequ + where takePadded sublen subsequ = case take sublen subsequ of + subsubsequ | length subsubsequ == sublen -> subsubsequ + | otherwise -> subsubsequ ++ replicate (sublen - length subsubsequ) 'N' + start = fromIntegral istart + len = fromIntegral ilen + +instance SeqLike BS.ByteString where + length = fromIntegral . BS.length + ntAt sequ ipos | pos >= 0 && pos < BS.length sequ = Just $ BS.index sequ pos + | otherwise = Nothing + where pos = fromIntegral ipos + subseq = bsSubseq + subseqPad = bsSubseqPad + concat = BS.concat + +bsSubseq :: (Integral n, Integral m) => n -> m -> BS.ByteString -> Maybe BS.ByteString +bsSubseq istart ilen sequ + | start < 0 || start + len > BS.length sequ = Nothing + | otherwise = Just . BS.take len . BS.drop start $ sequ + where start = fromIntegral istart + len = fromIntegral ilen + +bsSubseqPad :: (Integral n, Integral m) => n -> m -> BS.ByteString -> BS.ByteString +bsSubseqPad istart ilen sequ + | start + len <= 0 = BS.replicate len 'N' + | start >= BS.length sequ = BS.replicate len 'N' + | start < 0 = BS.replicate (negate start) 'N' `BS.append` takePadded (len + start) sequ + | otherwise = takePadded len $ BS.drop start sequ + where takePadded sublen subsequ + | sublen <= BS.length subsequ = BS.take sublen subsequ + | otherwise = subsequ `BS.append` BS.replicate (sublen - BS.length subsequ) 'N' + start = fromIntegral istart + len = fromIntegral ilen + +instance SeqLike LBS.ByteString where + length = fromIntegral . LBS.length + ntAt = lbsNtAt + subseq = lbsSubseq + subseqPad = lbsSubseqPad + concat = LBS.concat + +lbsNtAt :: (Integral n) => LBS.ByteString -> n -> Maybe Char +lbsNtAt sequ pos | pos < 0 = Nothing + | otherwise = case LBS.drop (fromIntegral pos) sequ of + trimmed | LBS.null trimmed -> Nothing + | otherwise -> Just . LBS.head $ trimmed + +lbsSubseq :: (Integral n, Integral m) => n -> m -> LBS.ByteString -> Maybe LBS.ByteString +lbsSubseq istart ilen sequ + | start < 0 = Nothing + | otherwise = case LBS.take len . LBS.drop start $ sequ of + subsequ | LBS.length subsequ < len -> Nothing + | otherwise -> Just subsequ + where start = fromIntegral istart + len = fromIntegral ilen + +lbsSubseqPad :: (Integral n, Integral m) => n -> m -> LBS.ByteString -> LBS.ByteString +lbsSubseqPad istart ilen sequ + | start + len <= 0 = LBS.replicate len 'N' + | start < 0 = LBS.replicate (negate start) 'N' `LBS.append` takePadded (len + start) sequ + | otherwise = takePadded len $ LBS.drop start sequ + where takePadded sublen subsequ = case LBS.take sublen subsequ of + subsubsequ | LBS.length subsubsequ == sublen -> subsubsequ + | otherwise -> subsubsequ `LBS.append` LBS.replicate (sublen - LBS.length subsubsequ) 'N' + start = fromIntegral istart + len = fromIntegral ilen diff --git a/src/Bio/SeqLoc/SpliceLocation.hs b/src/Bio/SeqLoc/SpliceLocation.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/SpliceLocation.hs @@ -0,0 +1,188 @@ +{-# LANGUAGE TypeFamilies, FlexibleContexts, OverloadedStrings #-} + +{-| Data type for a more general sequence location consiting of +disjoint ranges of positions on a sequence. + +Throughout, /sequence position/ refers to a 'Pos.Pos' which includes a +strand, as opposed to an /offset/, which refers to a 'Pos.Offset' with +no strand. + -} + +module Bio.SeqLoc.SpliceLocation ( + -- * Sequence locations + SpliceLoc + , fromContigs + , locOutof, locInto + , mergeContigs, mergeAdjContigs + ) where + +import Prelude hiding (length) + +import Control.Applicative +import Control.Arrow ((***)) +import Control.Monad +import qualified Data.ByteString.Char8 as BS +import Data.List (foldl') +import Data.Maybe + +import qualified Data.Attoparsec.Zepto as ZP + +import Bio.SeqLoc.LocRepr +import Bio.SeqLoc.Location +import qualified Bio.SeqLoc.Position as Pos +import Bio.SeqLoc.Strand +import qualified Bio.SeqLoc.SeqLike as SeqLike + +-- | General (disjoint) sequence region consisting of a concatenated +-- set of one or more contiguous regions. +data SpliceLoc = SpliceLocLast { contig :: !ContigLoc } + | SpliceLocPrev { contig :: !ContigLoc + , _next :: !SpliceLoc + } + deriving (Eq, Ord, Show) + +singleton :: ContigLoc -> SpliceLoc +singleton = SpliceLocLast + +consContig :: ContigLoc -> SpliceLoc -> Maybe SpliceLoc +consContig c sploc | goodJunction c (contig sploc) = Just $! SpliceLocPrev c sploc + | otherwise = Nothing + +goodJunction :: ContigLoc -> ContigLoc -> Bool +goodJunction c5 c3 = sameStrand && inOrder + where sameStrand = strand c5 == strand c3 + c5end = snd . bounds $ c5 + c3start = fst . bounds $ c3 + inOrder = case strand c5 of + Fwd -> c5end < c3start + RevCompl -> c5end > c3start + +fromContigs :: [ContigLoc] -> Maybe SpliceLoc +fromContigs [] = Nothing +fromContigs [l] = Just $! singleton l +fromContigs (l:rest) = fromContigs rest >>= consContig l + +tails :: SpliceLoc -> [SpliceLoc] +tails sll@(SpliceLocLast _) = [sll] +tails slp@(SpliceLocPrev _ n) = slp : (tails n) + +instance Stranded SpliceLoc where + revCompl sll = fromMaybe badRevCompl . foldl' addprev newlast $ slrest + where (sl0:slrest) = tails sll + newlast = Just $! singleton . revCompl . contig $ sl0 + addprev slnew slold = slnew >>= consContig (revCompl . contig $ slold) + badRevCompl = error $ "Bad junction doing reverse complement on " ++ (BS.unpack . repr) sll + +instance LocRepr SpliceLoc where + repr = BS.intercalate (BS.singleton ';') . map repr . contigs + unrepr = (fromContigs <$> scan) >>= maybe (fail "bad contig order") return + where scan = liftA2 (:) unrepr ((ZP.string ";" *> scan) <|> pure []) + +instance Location SpliceLoc where + strand = strand . contig + length = foldl' (\len c -> len + length c) 0 . contigs + bounds = (minimum *** maximum) . unzip . map bounds . contigs + seqData sequ = liftM SeqLike.concat . mapM (seqData sequ) . contigs + seqDataPad sequ = SeqLike.concat . map (seqDataPad sequ) . contigs + posInto pos = posIntoContigs pos . contigs + posOutof pos = posOutofContigs pos . contigs + startPos = startPos . firstContig + endPos = endPos . lastContig + clocInto = slocClocInto + clocOutof = slocClocOutof + extend = slocExtend + offsetWithin off = or . map (offsetWithin off) . contigs + posWithin pos = or . map (posWithin pos) . contigs + contigOverlaps c = any (contigOverlaps c ) . contigs + toContigs = contigs + +locOutof :: (Location l) => SpliceLoc -> l -> Maybe SpliceLoc +locOutof sploc outer = mapM (flip clocOutof outer) (toContigs sploc) >>= + fromContigs . concat . map toContigs + +locInto :: (Location l) => SpliceLoc -> l -> Maybe SpliceLoc +locInto sploc outer = mapM (flip clocInto outer) (toContigs sploc) >>= + fromContigs . mergeContigs . concat . map toContigs + +mergeContigs :: [ContigLoc] -> [ContigLoc] +mergeContigs [] = [] +mergeContigs [clast] = [clast] +mergeContigs (cprev:rest@(_:_)) = case mergeContigs rest of + [] -> error $ "mergeContigs: empty rest" + mergerest@(cnext:afternext) -> case mergeAdjContigs cprev cnext of + Just cmerge -> cmerge : afternext + Nothing -> cprev : mergerest + +mergeAdjContigs :: ContigLoc -> ContigLoc -> Maybe ContigLoc +mergeAdjContigs clocprev clocnext + | startPos clocnext == endPos (extend (0, 1) clocprev) + = Just $! fromPosLen (startPos clocprev) (length clocprev + length clocnext) + | otherwise = Nothing + +contigs :: SpliceLoc -> [ContigLoc] +contigs (SpliceLocLast c) = [c] +contigs (SpliceLocPrev c n) = c : contigs n + +firstContig :: SpliceLoc -> ContigLoc +firstContig = contig + +lastContig :: SpliceLoc -> ContigLoc +lastContig (SpliceLocLast c) = c +lastContig (SpliceLocPrev _ n) = lastContig n + +--reloffsets :: SpliceLoc -> [Pos.Offset] +--reloffsets = init . scanl (\off c -> off + length c) 0 . contigs + +contigsAndOffsets :: SpliceLoc -> [(ContigLoc, Pos.Offset)] +contigsAndOffsets = go 0 + where go off (SpliceLocLast c) = [(c, off)] + go off (SpliceLocPrev c n) = (c, off) : (go (off + length c) n) + +posIntoContigs :: Pos.Pos -> [ContigLoc] -> Maybe Pos.Pos +posIntoContigs pos = go 0 + where go _ [] = Nothing + go dlen (c:rest) = maybe onRest onContig . posInto pos $ c + where onContig pin = Just $! pin `Pos.slide` dlen + onRest = go (dlen + length c) rest + +posOutofContigs :: Pos.Pos -> [ContigLoc] -> Maybe Pos.Pos +posOutofContigs _ [] = Nothing +posOutofContigs pos (c:rest) = maybe onRest onContig . posOutof pos $ c + where onRest = posOutofContigs (Pos.slide pos . negate . length $ c) rest + onContig = Just + +slocClocInto :: ContigLoc -> SpliceLoc -> Maybe ContigLoc +slocClocInto subcloc = listToMaybe . catMaybes . map into . contigsAndOffsets + where into (cloc, off) = liftM (slide off) . clocInto subcloc $ cloc + +slocClocOutof :: ContigLoc -> SpliceLoc -> Maybe SpliceLoc +slocClocOutof cloc + = liftM (stranded (strand cloc) . fromContigsErr) . regionOutof (offset5 cloc) (length cloc) . contigs + where fromContigsErr ctgs = fromMaybe badContigs $! fromContigs ctgs + where badContigs = error . unwords $ + [ "bad contig order in slocClocOutof" ] ++ + map (BS.unpack . repr) ctgs + + +regionOutof :: Pos.Offset -> Pos.Offset -> [ContigLoc] -> Maybe [ContigLoc] +regionOutof _ _ [] = Nothing +regionOutof pos5 len (cloc0:rest) + | pos5 < 0 = Nothing + | len < 0 = error $ "Bio.SeqLoc.SpliceLocation.regionOutof: input, " ++ show (pos5, len, cloc0) + | pos5 >= len0 = regionOutof (pos5 - len0) len rest + | (pos5 + len <= len0) = let subcloc = fromPosLen (Pos.Pos pos5 Fwd) len + in case clocOutof subcloc cloc0 of + Just out0 -> Just [out0] + Nothing -> error $ "regionOutof: final bounds failure, " ++ (BS.unpack . repr) subcloc + | otherwise = let outlen0 = len0 - pos5 + subcloc = fromPosLen (Pos.Pos pos5 Fwd) outlen0 + in case clocOutof subcloc cloc0 of + Just out0 -> liftM (out0 :) . regionOutof 0 (len - outlen0) $ rest + Nothing -> error $ "regionOutof: internal bounds failure, " ++ (BS.unpack . repr) subcloc + where len0 = length cloc0 + +slocExtend :: (Pos.Offset, Pos.Offset) -> SpliceLoc -> SpliceLoc +slocExtend (ext5, ext3) sloc = extended3 { contig = extend (ext5, 0) . contig $ extended3 } + where extend3 (SpliceLocPrev c n) = SpliceLocPrev c (extend3 n) + extend3 (SpliceLocLast c) = SpliceLocLast $ extend (0, ext3) c + extended3 = extend3 sloc diff --git a/src/Bio/SeqLoc/Strand.hs b/src/Bio/SeqLoc/Strand.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/Strand.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE OverloadedStrings, TypeSynonymInstances #-} + +{-| Utilities for manipulating nucleotide sequences and locations on +nucleotide sequences that occur on a forward or a reverse-complement +strand. + +-} + +module Bio.SeqLoc.Strand ( Strand(..) + , compl + , Stranded(..), stranded + ) + where + +import Control.Applicative +import Data.ByteString.Internal (c2w, w2c) +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as LBS +import Data.Word (Word8) + +import qualified Data.Attoparsec.Zepto as ZP + +import Bio.SeqLoc.LocRepr + +-- | Complement of a nucleotide character, swap A/T and G/C preserving +-- case and leave all other characters unchanged. +compl :: Char -> Char +compl 'a' = 't' +compl 'c' = 'g' +compl 'g' = 'c' +compl 't' = 'a' +compl 'A' = 'T' +compl 'C' = 'G' +compl 'G' = 'C' +compl 'T' = 'A' +compl ch = ch + +-- | Sequence strand +data Strand = Fwd | RevCompl deriving (Eq, Ord, Show, Read, Bounded, Enum) + +instance LocRepr Strand where + repr Fwd = BS.pack "(+)" + repr RevCompl = BS.pack "(-)" + unrepr = (ZP.string "(+)" *> return Fwd) <|> + (ZP.string "(-)" *> return RevCompl) + +-- | A nucleotide sequence or location on a nucleotide sequence that +-- lies on a specific strand and has an orientation. +class Stranded s where + revCompl :: s -> s + +-- | Convert the orientation of a 'Stranded' thing based on a +-- specified 'Strand' +stranded :: (Stranded s) => Strand -> s -> s +stranded Fwd = id +stranded RevCompl = revCompl + +instance Stranded Strand where + revCompl Fwd = RevCompl + revCompl RevCompl = Fwd + +instance Stranded Char where + revCompl = compl + +instance Stranded Word8 where + revCompl = c2w . compl . w2c + +instance Stranded String where + revCompl = reverse . map compl + +instance Stranded LBS.ByteString where + revCompl = LBS.reverse . LBS.map compl + +instance Stranded BS.ByteString where + revCompl = BS.reverse . BS.map compl + diff --git a/src/Bio/SeqLoc/Transcript.hs b/src/Bio/SeqLoc/Transcript.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/Transcript.hs @@ -0,0 +1,116 @@ +module Bio.SeqLoc.Transcript + ( + -- * Type for splice junctions + Junction (..) + , fromDonorAcceptor, donor, acceptor + , junctions + -- * Representation of transcript + , Transcript(..), utr5, utr3 + , cdsLocation + , sortContigs + ) + where + +import Control.Applicative +import Control.Monad +import qualified Data.ByteString.Char8 as BS +import Data.List +import Data.Ord + +import qualified Data.Attoparsec.Zepto as ZP + +import Bio.SeqLoc.LocRepr +import qualified Bio.SeqLoc.Location as Loc +import Bio.SeqLoc.OnSeq +import qualified Bio.SeqLoc.Position as Pos +import qualified Bio.SeqLoc.SpliceLocation as SpLoc +import Bio.SeqLoc.Strand + +-- | Splice junctions, which are isomorphic to the introns they span, +-- but which support other biologically relevant constructors and +-- accessors. +newtype Junction = Junction { intron :: Loc.ContigLoc } deriving (Show) + +slash :: BS.ByteString +slash = BS.pack "/" + +instance LocRepr Junction where + repr j = BS.concat [ repr . donor $ j, slash, repr . acceptor $ j ] + unrepr = fromDonorAcceptor <$> unrepr <*> (ZP.string slash *> unrepr) + +-- | Create a splice junction from a donor position (the last position +-- in the 5' exon) and the acceptor position (the first position in +-- the 3' exon). +fromDonorAcceptor :: Pos.Pos -> Pos.Pos -> Junction +fromDonorAcceptor d a = let len = 1 + abs (Pos.offset a - Pos.offset d) + in case Pos.strand d of + Fwd -> Junction $! Loc.fromPosLen (Pos.slide d 1) len + RevCompl -> Junction $! Loc.fromPosLen (Pos.slide d (-1)) len + +-- | Donor position, i.e., the last position in the 5' exon around a +-- junction. +donor :: Junction -> Pos.Pos +donor = Loc.startPos . Loc.extend (1, 0) . intron + +-- | Acceptor position, i.e., the first position in the 3' exon around +-- a junction. +acceptor :: Junction -> Pos.Pos +acceptor = Loc.endPos . Loc.extend (0, 1) . intron + +-- | List of splice junctions from a spliced location, in order. +junctions :: SpLoc.SpliceLoc -> [Junction] +junctions sploc = zipWith junction contigs (drop 1 contigs) + where contigs = Loc.toContigs sploc + junction c5 c3 = let p5 = Loc.endPos . Loc.extend (0, 1) $ c5 + p3 = Loc.startPos . Loc.extend (1, 0) $ c3 + len = 1 + abs (Pos.offset p3 - Pos.offset p5) + in Junction $ Loc.fromPosLen p5 len + + + +-- | Representation of a genomic transcript, with a gene and a +-- transcript identifier, along with the genomic location of the +-- processed transcript and an optional coding sequence on that +-- transcript. +data Transcript = Transcript { geneId :: !SeqName -- ^ Gene or locus name for a collection of transcripts + , trxId :: !SeqName -- ^ Specific transcript identifier + , location :: !SpliceSeqLoc -- ^ Sequence location of processed transcript + , cds :: !(Maybe Loc.ContigLoc) -- ^ Location of CDS on the transcript + } deriving (Show) + +-- | 'Just' the location of the 5' UTR on the transcript, or 'Nothing' +-- if there is no 'cds' on the transcript or if the 'cds' location +-- begins at the first nucleotide of the transcript--if a region is +-- returned it will have positive length. +utr5 :: Transcript -> Maybe Loc.ContigLoc +utr5 trx = cds trx >>= utr5loc + where utr5loc cdsloc = case Loc.startPos cdsloc of + (Pos.Pos startoff Fwd) | startoff > 0 -> Just $! Loc.fromBoundsStrand 0 (startoff - 1) Fwd + _ -> Nothing + +-- | 'Just' the location of the 3' UTR on the transcript, or 'Nothing' +-- if there is no 'cds' on the transcript or if the 'cds' location +-- ends at the last nucleotide of the transcript--if a region is +-- returned it will have positive length. +utr3 :: Transcript -> Maybe Loc.ContigLoc +utr3 trx = cds trx >>= utr3loc + where utr3loc cdsloc = case Loc.endPos cdsloc of + (Pos.Pos endoff Fwd) | endoff < trxlast -> Just $! Loc.fromBoundsStrand (endoff + 1) trxlast Fwd + _ -> Nothing + trxlast = snd . Loc.bounds . unOnSeq . location $ trx + +-- | Genomic location of CDS within the transcript +cdsLocation :: Transcript -> Maybe SpliceSeqLoc +cdsLocation trx = cds trx >>= liftM (OnSeq name) . flip Loc.clocOutof loc + where (OnSeq name loc) = location trx + +-- | 'Just' the input contigs sorted in stranded order, when all lie +-- on the same strand, or 'Nothing' if they are not all on the same +-- strand. +sortContigs :: [Loc.ContigLoc] -> Maybe [Loc.ContigLoc] +sortContigs [] = Nothing +sortContigs cs@(c0:_)= liftM sortStrand contigStrand + where contigStrand | all ((== Loc.strand c0) . Loc.strand) cs = Just . Loc.strand $ c0 + | otherwise = Nothing + sortStrand Fwd = sortBy (comparing Loc.offset5) cs + sortStrand RevCompl = sortBy (comparing (negate . Loc.offset5)) cs diff --git a/test/TestMain.hs b/test/TestMain.hs new file mode 100644 --- /dev/null +++ b/test/TestMain.hs @@ -0,0 +1,392 @@ +{-# LANGUAGE ExistentialQuantification #-} +module Main + where + +import Control.Monad +import qualified Data.ByteString as BSW +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as LBS +import Data.ByteString.Internal (c2w, w2c) +import Data.Char +import Data.Either +import Data.Ix (inRange) +import Data.List +import Data.Maybe +import System.Random + +import Test.QuickCheck + +import Bio.SeqLoc.LocRepr +import qualified Bio.SeqLoc.Location as Loc +import Bio.SeqLoc.OnSeq +import qualified Bio.SeqLoc.Position as Pos +import qualified Bio.SeqLoc.SpliceLocation as SpLoc +import Bio.SeqLoc.Strand +import qualified Bio.SeqLoc.SeqLike as SeqLike + +main :: IO () +main = mapM_ runTest tests + +tests :: [Test] +tests = [ T "Strand revCompl" test_Strand_revCompl + , T "Char revCompl" property_Char_revCompl + , T "ByteString revCompl" property_ByteString_revCompl + , T "Sequence revCompl" property_Sequence_revCompl + + , T "Pos revCompl" test_Pos_revCompl + , T "Pos atPos" property_Pos_atPos + , T "Pos atPos2" property_Pos_atPos2 + , T "Pos repr" test_Pos_repr + + , T "Contig revCompl" test_Contig_RevCompl + , T "Contig pos into/outof inverse" property_ContigIntoOutof + , T "Contig pos outof/into inverse" property_ContigOutofInto + , T "Contig loc into/outof inverse" property_ContigLocIntoOutof + , T "Contig loc outof/into inverse" property_ContigLocOutofInto + , T "Contig into based on bounds" test_Contig_IntoBounds + , T "Contig outof based on bounds" test_Contig_OutofBounds + , T "Contig seqData" property_Contig_seqData + , T "Contig seqDataPadded" property_Contig_seqDataPadded + , T "Contig seqData2" property_Contig_seqData2 + , T "Contig extend/revCompl" property_Contig_extendRevCompl + , T "Contig fromStartEnd" property_Contig_fromStartEnd + , T "Contig fromBoundsStrand" property_Contig_fromBoundsStrand + , T "Contig overlaps" property_Contig_overlaps + , T "Contig repr" test_Contig_repr + + , T "Loc revCompl" test_Loc_RevCompl + , T "Loc pos into/outof inverse" property_LocIntoOutof + , T "Loc pos outof/into inverse" property_LocOutofInto + , T "Loc outof based on bounds" test_Loc_OutofBounds + , T "Loc loc outof/into inverse" property_LocCLocOutofInto + , T "Loc outof association" property_LocOutofAssoc + , T "Loc locOutof" property_SpLocOutof + , T "Loc locOutof valid" property_SpLocOutofGood + , T "Loc within" property_Loc_Within + , T "Loc seqData" property_Loc_seqData + , T "Loc seqDataPadded" property_Loc_seqDataPadded + , T "SpLoc seqData2" property_SpLoc_seqData2 + , T "SpLoc repr" test_SpLoc_repr + , T "SpLoc termini/revCompl" property_SpLoc_terminiRevCompl + , T "SpLoc termini/extend" property_SpLoc_terminiExtend + ] + + +-- Bio.BioSeq.Stranded + +genNtByteString :: Int -> Gen BS.ByteString +genNtByteString = liftM BS.pack . flip replicateM (elements "ACGT") + +genName :: Gen SeqName +genName = liftM (SeqName . BS.pack) $ genNameLength >>= flip replicateM genNameChar + where genNameLength = choose (1, 20) + genNameChar = elements $ ['a'..'z'] ++ ['A'..'Z'] ++ ['0'..'9'] ++ "-_" + +instance Arbitrary SeqName where + arbitrary = genName + +test_revCompl :: (Eq s, Stranded s) => s -> Bool +test_revCompl s = (revCompl . revCompl) s == s + +test_repr :: (LocRepr l, Eq l) => l -> Bool +test_repr l = (unreprMaybe . repr $ l) == Just l + +test_Strand_revCompl :: Strand -> Bool +test_Strand_revCompl = test_revCompl + +property_Char_revCompl :: Property +property_Char_revCompl = forAll (elements "ACGTacgtnN") test_revCompl + +property_ByteString_revCompl :: Property +property_ByteString_revCompl = forAll (sized genNtByteString) test_revCompl + +property_Sequence_revCompl :: Property +property_Sequence_revCompl + = forAll arbitrary $ \name -> + let mkSeq = OnSeq name + in forAll (sized genNtByteString) $ \sequ -> + (unOnSeq . revCompl . mkSeq) sequ == revCompl sequ + +-- Bio.BioSeq.Position + +test_Pos_revCompl :: Pos.Pos -> Bool +test_Pos_revCompl = test_revCompl + +property_Pos_atPos :: Pos.Pos -> Property +property_Pos_atPos pos + = forAll genPositiveOffset $ \seqlen -> + forAll (genNtByteString $ fromIntegral seqlen) $ \sequ -> + let actual = Pos.atPos sequ pos + in if and [ Pos.offset pos >= 0, Pos.offset pos < seqlen ] + then let fwdNt = BS.index sequ (fromIntegral . Pos.offset $ pos) + in case Pos.strand pos of + Fwd -> actual == Just fwdNt + RevCompl -> actual == Just (compl $ fwdNt) + else actual == Nothing + +property_Pos_atPos2 :: Pos.Pos -> Property +property_Pos_atPos2 pos + = forAll genPositiveOffset $ \seqlen -> + forAll (genNtByteString $ fromIntegral seqlen) $ \sequ -> + and [ Pos.atPos sequ pos == Pos.atPos (LBS.fromChunks [sequ]) pos + , Pos.atPos sequ pos == Pos.atPos (BS.unpack sequ) pos + ] + +test_Pos_repr :: Pos.Pos -> Bool +test_Pos_repr = test_repr + +-- Bio.BioSeq.Location + +instance Arbitrary Strand where + arbitrary = elements [Fwd, RevCompl] + +instance Arbitrary Pos.Pos where + arbitrary = liftM2 Pos.Pos genOffset arbitrary + +instance Arbitrary Loc.ContigLoc where + arbitrary = liftM2 Loc.fromPosLen arbitrary genPositiveOffset + +test_Contig_RevCompl :: Loc.ContigLoc -> Bool +test_Contig_RevCompl = test_revCompl + +property_ContigIntoOutof :: Loc.ContigLoc -> Pos.Pos -> Property +property_ContigIntoOutof contig pos + = let !mInpos = Loc.posInto pos contig + !mOutpos = mInpos >>= flip Loc.posOutof contig + in (isJust mInpos) ==> mOutpos == (Just pos) + +property_ContigOutofInto :: Loc.ContigLoc -> Pos.Pos -> Property +property_ContigOutofInto contig pos + = let !mOutpos = Loc.posOutof pos contig + !mInpos = mOutpos >>= flip Loc.posInto contig + in (isJust mOutpos) ==> mInpos == (Just pos) + +property_ContigLocIntoOutof :: Loc.ContigLoc -> Loc.ContigLoc -> Property +property_ContigLocIntoOutof subcloc supercloc + = let !mIncloc = Loc.clocInto subcloc supercloc + !mOutcloc = mIncloc >>= flip Loc.clocOutof supercloc + in (isJust mIncloc) ==> mOutcloc == (Just subcloc) + +property_ContigLocOutofInto :: Loc.ContigLoc -> Loc.ContigLoc -> Property +property_ContigLocOutofInto subcloc supercloc + = let !mOutcloc = Loc.clocOutof subcloc supercloc + !mIncloc = mOutcloc >>= flip Loc.clocInto supercloc + in (isJust mOutcloc) ==> mIncloc == (Just subcloc) + +test_Contig_IntoBounds :: Loc.ContigLoc -> Pos.Pos -> Bool +test_Contig_IntoBounds contig pos + = let !mInpos = Loc.posInto pos contig + !offset = Pos.offset pos + !(cstart, cend) = Loc.bounds contig + in (isJust mInpos) == (offset >= cstart && offset <= cend) + +test_Contig_OutofBounds :: Loc.ContigLoc -> Pos.Pos -> Bool +test_Contig_OutofBounds contig pos + = let !offset = Pos.offset pos + in (isJust $ Loc.posOutof pos contig) == (offset >= 0 && offset < Loc.length contig) + +property_Contig_seqData :: Loc.ContigLoc -> Property +property_Contig_seqData contig + = forAll (genNonNegOffset >>= genNtByteString . fromIntegral) $ \sequ -> + let seqData = Loc.seqData sequ contig + padded = Loc.seqDataPad sequ contig + in case seqData of + (Just subsequ) -> and [ padded == subsequ, 'N' `BS.notElem` padded ] + Nothing -> 'N' `BS.elem` padded + +property_Contig_seqDataPadded :: Loc.ContigLoc -> Property +property_Contig_seqDataPadded contig + = forAll (genNonNegOffset >>= genNtByteString . fromIntegral) $ \sequ -> + (BS.pack $ map (fromMaybe 'N' . Pos.atPos sequ) contigPoses) == Loc.seqDataPad sequ contig + where contigPoses = mapMaybe (flip Loc.posOutof contig . flip Pos.Pos Fwd) [0..(Loc.length contig - 1)] + +property_Contig_seqData2 :: Loc.ContigLoc -> Property +property_Contig_seqData2 contig + = forAll (genNonNegOffset >>= genNtByteString . fromIntegral) $ \sequ -> + let toLBS = LBS.fromChunks . (: []) + fromLBS = BS.concat . LBS.toChunks + in and [ Loc.seqData sequ contig == liftM fromLBS (Loc.seqData (toLBS sequ) contig) + , Loc.seqDataPad sequ contig == fromLBS (Loc.seqDataPad (toLBS sequ) contig) + , Loc.seqData sequ contig == liftM BS.pack (Loc.seqData (BS.unpack sequ) contig) + , Loc.seqDataPad sequ contig == BS.pack (Loc.seqDataPad (BS.unpack sequ) contig) + ] + +property_Contig_extendRevCompl :: Loc.ContigLoc -> Property +property_Contig_extendRevCompl contig + = forAll (liftM2 (,) genNonNegOffset genNonNegOffset) $ \(ext5, ext3) -> + (revCompl $ Loc.extend (ext5, ext3) contig) == (Loc.extend (ext3, ext5) $ revCompl contig) + +property_Contig_fromStartEnd :: Loc.ContigLoc -> Property +property_Contig_fromStartEnd contig + = (Loc.length contig > 1) ==> + (Loc.fromStartEnd (Pos.offset $ Loc.startPos contig) (Pos.offset $ Loc.endPos contig)) == contig + +property_Contig_fromBoundsStrand :: Loc.ContigLoc -> Property +property_Contig_fromBoundsStrand contig + = (Loc.length contig > 1) ==> + (Loc.fromBoundsStrand (fst . Loc.bounds $ contig) (snd . Loc.bounds $ contig) (Loc.strand contig)) == contig + +property_Contig_overlaps :: Loc.ContigLoc -> Loc.ContigLoc -> Bool +property_Contig_overlaps cloc1 cloc2 + = (cloc1 `Loc.contigOverlaps` cloc2) == + and [ Loc.strand cloc1 == Loc.strand cloc2 + , or [ isJust . Loc.posInto (Loc.startPos cloc1) $ cloc2 + , isJust . Loc.posInto (Loc.startPos cloc2) $ cloc1 + , isJust . Loc.posInto (Loc.endPos cloc1) $ cloc2 + , isJust . Loc.posInto (Loc.endPos cloc2) $ cloc1 + ] + ] + +test_Contig_repr :: Loc.ContigLoc -> Bool +test_Contig_repr = test_repr + +-- Bio.BioSeq.Location + +genInvertibleLoc :: Gen SpLoc.SpliceLoc +genInvertibleLoc = sized $ \sz -> do ncontigs <- choose (1, sz + 1) + fwdloc <- liftM (fromJust . SpLoc.fromContigs) + $ genContigs ncontigs + rc <- arbitrary + if rc then return $ revCompl fwdloc else return fwdloc + where genContigs = liftM (reverse . foldl' intervalsToContigs []) . genIntervals + genIntervals nints = replicateM nints $ liftM2 (,) genPositiveOffset genPositiveOffset + intervalsToContigs [] (init5, len) = [Loc.fromPosLen (Pos.Pos init5 Fwd) len] + intervalsToContigs prevs@(prev:_) (nextoffset, nextlen) + = let !prevend = Loc.offset5 prev + Loc.length prev + in (Loc.fromPosLen (Pos.Pos (prevend + nextoffset) Fwd) nextlen):prevs + +instance Arbitrary SpLoc.SpliceLoc where + arbitrary = genInvertibleLoc + +test_Loc_RevCompl :: SpLoc.SpliceLoc -> Bool +test_Loc_RevCompl = test_revCompl + +property_LocIntoOutof :: SpLoc.SpliceLoc -> Pos.Pos -> Property +property_LocIntoOutof loc pos + = let !mInpos = Loc.posInto pos loc + !mOutpos = mInpos >>= flip Loc.posOutof loc + in (isJust mInpos) ==> mOutpos == (Just pos) + +property_LocOutofInto :: Pos.Pos -> Property +property_LocOutofInto pos + = forAll genInvertibleLoc $ \loc -> + let !mOutpos = Loc.posOutof pos loc + !mInpos = mOutpos >>= flip Loc.posInto loc + in (isJust mOutpos) ==> mInpos == (Just pos) + +test_Loc_OutofBounds :: SpLoc.SpliceLoc -> Pos.Pos -> Bool +test_Loc_OutofBounds loc pos + = let !offset = Pos.offset pos + in (isJust $ Loc.posOutof pos loc) == (offset >= 0 && offset < Loc.length loc) + +property_LocCLocOutofInto :: Loc.ContigLoc -> Property +property_LocCLocOutofInto cloc + = forAll genInvertibleLoc $ \loc -> + let !mOutloc = Loc.clocOutof cloc loc + !mInloc = mOutloc >>= mapM (flip Loc.clocInto loc) . Loc.toContigs >>= return . fromJust . SpLoc.fromContigs + in (isJust mOutloc) ==> and [ liftM Loc.length mInloc == Just (Loc.length cloc) + , liftM Loc.bounds mInloc == Just (Loc.bounds cloc) + ] + +property_LocOutofAssoc :: SpLoc.SpliceLoc -> Loc.ContigLoc -> Pos.Pos -> Property +property_LocOutofAssoc loc cloc pos + = let !mOutloc = Loc.clocOutof cloc loc + !mOutpos = mOutloc >>= \outloc -> Loc.posOutof pos outloc + in (isJust mOutpos) ==> mOutpos == (Loc.posOutof pos cloc >>= \outpos -> Loc.posOutof outpos loc) + +property_SpLocOutof :: SpLoc.SpliceLoc -> SpLoc.SpliceLoc -> Bool +property_SpLocOutof subloc outerloc = + let !mOutofContigs = liftM Loc.toContigs $ SpLoc.locOutof subloc outerloc + !mContigsOutof = liftM (concat . map Loc.toContigs) $ + mapM (flip Loc.clocOutof outerloc) $ + Loc.toContigs subloc + in mOutofContigs == mContigsOutof + +property_SpLocOutofGood :: SpLoc.SpliceLoc -> SpLoc.SpliceLoc -> Property +property_SpLocOutofGood subloc outerloc = + let !mOutofContigs = liftM Loc.toContigs $ SpLoc.locOutof subloc outerloc + !mContigsOutof = liftM (concat . map Loc.toContigs) $ + mapM (flip Loc.clocOutof outerloc) $ + Loc.toContigs subloc + in (isJust mOutofContigs) ==> mOutofContigs == mContigsOutof + +property_Loc_seqData :: SpLoc.SpliceLoc -> Property +property_Loc_seqData loc + = forAll (genNonNegOffset >>= genNtByteString . fromIntegral) $ \sequ -> + let seqData = Loc.seqData sequ loc + padded = Loc.seqDataPad sequ loc + in case seqData of + (Just subsequ) -> and [ padded == subsequ, 'N' `BS.notElem` padded ] + Nothing -> 'N' `BS.elem` padded + +property_Loc_seqDataPadded :: SpLoc.SpliceLoc -> Property +property_Loc_seqDataPadded loc + = forAll (genNonNegOffset >>= genNtByteString . fromIntegral) $ \sequ -> + (BS.pack $ map (fromMaybe 'N' . Pos.atPos sequ) locPoses) == Loc.seqDataPad sequ loc + where locPoses = mapMaybe (flip Loc.posOutof loc . flip Pos.Pos Fwd) [0..(Loc.length loc - 1)] + +property_SpLoc_seqData2 :: SpLoc.SpliceLoc -> Property +property_SpLoc_seqData2 sploc + = forAll (genNonNegOffset >>= genNtByteString . fromIntegral) $ \sequ -> + let toLBS = LBS.fromChunks . (: []) + fromLBS = BS.concat . LBS.toChunks + in and [ Loc.seqData sequ sploc == liftM fromLBS (Loc.seqData (toLBS sequ) sploc) + , Loc.seqDataPad sequ sploc == fromLBS (Loc.seqDataPad (toLBS sequ) sploc) + , Loc.seqData sequ sploc == liftM BS.pack (Loc.seqData (BS.unpack sequ) sploc) + , Loc.seqDataPad sequ sploc == BS.pack (Loc.seqDataPad (BS.unpack sequ) sploc) + ] + +property_Loc_Within :: Pos.Pos -> Property +property_Loc_Within pos + = forAll genInvertibleLoc $ \loc -> + and [ (pos `Loc.posWithin` loc) == (maybe False ((/= RevCompl) . Pos.strand) $ Loc.posInto pos loc) + , ((Pos.offset pos) `Loc.offsetWithin` loc) == (isJust . Loc.posInto pos $ loc) + ] + +test_SpLoc_repr :: SpLoc.SpliceLoc -> Bool +test_SpLoc_repr = test_repr + +property_SpLoc_terminiRevCompl :: SpLoc.SpliceLoc -> Bool +property_SpLoc_terminiRevCompl loc + = and [ revCompl (Loc.startPos loc) == Loc.endPos (revCompl loc) + , revCompl (Loc.endPos loc) == Loc.startPos (revCompl loc) + ] + +property_SpLoc_terminiExtend :: SpLoc.SpliceLoc -> Property +property_SpLoc_terminiExtend loc + = forAll genNonNegOffset $ \ext5 -> + forAll genNonNegOffset $ \ext3 -> + let extloc = Loc.extend (ext5, ext3) loc + strandSlide pos doff = case Pos.strand pos of + Fwd -> Pos.slide pos doff + RevCompl -> Pos.slide pos (negate doff) + in and [ Loc.startPos extloc == strandSlide (Loc.startPos loc) (negate ext5) + , Loc.endPos extloc == strandSlide (Loc.endPos loc) ext3 + ] + +-- Utilities + +data Test = forall t . Testable t => T String t + +runTest :: Test -> IO () +runTest (T name test) = do + putStr $ name ++ replicate (40 - length name) '.' ++ " " + quickCheckWith args test + where args = stdArgs { maxDiscard = 100000 } + +-- | Constrained position generators + +genOffset :: Gen Pos.Offset +genOffset = do isneg <- arbitrary + nnoff <- genNonNegOffset + return $ (if isneg then negate else id) nnoff + +genNonNegOffset :: Gen Pos.Offset +genNonNegOffset = liftM (subtract 1) genPositiveOffset + +genPositiveOffset :: Gen Pos.Offset +genPositiveOffset = do scale <- chooseInteger (1, 10) + liftM fromIntegral $ chooseInteger (1, 2^scale) + where chooseInteger :: (Integer, Integer) -> Gen Integer + chooseInteger = choose +