packages feed

BiobaseTypes-0.2.1.0: Biobase/Types/Location.hs

-- | Annotate the genomic @Location@ of features or elements. A @Location@ is
-- always contiguous, using strand, 0-based position, and length.
-- Transformation to different systems of annotation is made possible.

module Biobase.Types.Location where

import Control.DeepSeq
import Control.Lens hiding (Index, index)
import Data.Coerce
import Data.Data
import Data.Data.Lens
import GHC.Generics (Generic)
import GHC.TypeNats
import Prelude hiding (length)
import qualified Data.ByteString as BS
import qualified Streaming.Internal as SI
import qualified Streaming.Prelude as SP
import Text.Printf

import Biobase.Types.BioSequence
import Biobase.Types.Index
import Biobase.Types.Position
import Biobase.Types.Strand
import Data.Info




-- | Operations on locations.

class ModifyLocation posTy seqTy where
  -- | Append to the left.
  locAppendLeft  :: seqTy -> Location i posTy seqTy -> Location i posTy seqTy
  -- | Append to the right.
  locAppendRight :: seqTy -> Location i posTy seqTy -> Location i posTy seqTy
  -- | Split a location.
  locSplitAt  :: Int -> Location i posTy seqTy -> (Location i posTy seqTy, Location i posTy seqTy)
  -- | Length of location
  locLength :: Location i posTy seqTy -> Int

locTake k = fst . locSplitAt k

locTakeEnd k loc = let l = locLength loc in snd $ locSplitAt (l-k) loc

locDrop k = snd . locSplitAt k

locDropEnd k loc = let l = locLength loc in fst $ locSplitAt (l-k) loc

locSplitEndAt k loc = let l = locLength loc in locSplitAt (l-k) loc



data Location ident posTy seqTy = Location
  { _locIdentifier  :: !(SequenceIdentifier ident)
  , _locPosition    :: !posTy
  , _locSequence    :: !seqTy
  }
  deriving stock (Show,Data,Typeable,Generic)
makeLenses ''Location

instance (NFData p, NFData s) => NFData (Location i p s)

retagLocation :: Location i posTy seqTy -> Location j posTy seqTy
{-# Inline retagLocation #-}
retagLocation = over locIdentifier coerce

instance ModifyLocation FwdPosition (BioSequence w) where
  {-# Inline locAppendLeft #-}
  locAppendLeft s loc = let l = s^._BioSequence.to BS.length
    in loc & locSequence %~ (s <>) & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (-. l) else p)
  {-# Inline locAppendRight #-}
  locAppendRight s loc = let l = s^._BioSequence.to BS.length
    in loc & locSequence %~ (<> s) & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~ (-. l) else p)
  {-# Inline locSplitAt #-}
  locSplitAt k loc =
    let (h',t') = loc^.locSequence._BioSequence.to (BS.splitAt k)
        hl = BS.length h' ; tl = BS.length t'
        h = loc & locSequence._BioSequence .~ h' & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~ (+. tl) else p)
        t = loc & locSequence._BioSequence .~ t' & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (+. hl) else p)
    in  (h,t)
  {-# Inline locLength #-}
  locLength = view (locSequence._BioSequence.to BS.length)

instance ModifyLocation FwdPosition Int where
  {-# Inline locAppendLeft #-}
  locAppendLeft k' loc = let k = max 0 $ min (loc^.locPosition.fwdStart.to toInt0) k' in
    loc & locSequence %~ (+ k) & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (-. k) else p)
  {-# Inline locAppendRight #-}
  locAppendRight k' loc = let k = max 0 $ min (loc^.locPosition.fwdStart.to toInt0) k' in
    loc & locSequence %~ (+ k) & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~ (-. k) else p)
  {-# Inline locSplitAt #-}
  locSplitAt k loc =
    let h' = max 0 . min k $ locLength loc
        t' = locLength loc - h'
        h = loc & locSequence .~ h' & locPosition %~ (\p -> if p^.fwdStrand == MinusStrand then p & fwdStart %~(+. t') else p)
        t = loc & locSequence .~ t' & locPosition %~ (\p -> if p^.fwdStrand == PlusStrand then p & fwdStart %~ (+. h') else p)
    in  (h,t)
  {-# Inline locLength #-}
  locLength = view locSequence

instance Reversing (Location i FwdPosition (BioSequence w)) where
  {-# Inline reversing #-}
  reversing = over (locSequence._BioSequence) BS.reverse . over (locPosition) reversing

instance Complement (BioSequence w) => Complement (Location i FwdPosition (BioSequence w)) where
  {-# Inline complement #-}
  complement = iso f f
    where f = over locSequence (view complement)

instance (Info (BioSequence w)) => Info (Location i FwdPosition (BioSequence w)) where
  info loc = printf "%s %s %s" (loc^.locIdentifier^.to show) (show $ loc^.locPosition) (loc^.locSequence.to info)

-- | Will extract a substring for a given biosequence. It is allowed to hand in partially or not at
-- all overlapping locational information. This will yield empty resulting locations.
--
-- This will convert the @FwdPosition@ strand, which in turn allows dealing with reverse-complement
-- searches.
--
-- @
-- 0123456789
--    3.3
-- @

subLocation :: Location i FwdPosition (BioSequence w) -> (FwdPosition, Int) -> Location i FwdPosition (BioSequence w)
{-# Inline subLocation #-}
subLocation s (p',l)
  | ss==PlusStrand = locTake l $ locDrop d s
  | ss==MinusStrand = locTakeEnd l $ locDropEnd d s
  where ss = s^.locPosition.fwdStrand
        p = if ss == p'^.fwdStrand then p' else reversing p'
        d = delta (s^.locPosition.fwdStart) (p^.fwdStart)

data PIS i p s = PIS
  { _pisPrefix  :: Maybe (Location i p s)
  , _pisInfix   :: !(Location i p s)
  , _pisSuffix  :: Maybe (Location i p s)
  }
  deriving stock (Show, Data)
makeLenses ''PIS

pis ifx = PIS Nothing ifx Nothing

retagPis :: PIS i p s -> PIS j p s
retagPis (PIS p i s) = PIS (fmap retagLocation p) (retagLocation i) (fmap retagLocation s)

-- | Given a @PIS@, this will return the @substring@ indicated by the location in the 2nd argument.
-- Allows for easy substring extraction, and retains the system of prefix/infix/suffix.
--
-- It is allowed to hand locations that only partially (or not at all) correspond to the @PIS@, but
-- then the resulting @PIS@ will be empty!

subPisLocation :: PIS i FwdPosition (BioSequence w) -> (FwdPosition, Int) -> PIS i FwdPosition (BioSequence w)
{-# Inline subPisLocation #-}
subPisLocation pis loc =
  let f z = subLocation z loc
  in  over (pisPrefix._Just) f . over pisInfix f $ over (pisSuffix._Just) f pis

instance (Reversing (Location i FwdPosition (BioSequence w))) => Reversing (PIS i FwdPosition (BioSequence w)) where
  {-# Inline reversing #-}
  reversing pis
    = over (pisPrefix._Just) reversing . over pisInfix reversing . over (pisSuffix._Just) reversing
    . set pisPrefix (pis^.pisSuffix) . set pisSuffix (pis^.pisPrefix) $ pis

instance Complement (BioSequence w) => Complement (PIS i FwdPosition (BioSequence w)) where
  {-# Inline complement #-}
  complement =
    let f = over pisInfix (view complement) . over (pisPrefix._Just) (view complement) . over (pisSuffix._Just) (view complement)
    in  iso f f

pisSequence :: Lens (PIS i p (BioSequence s)) (PIS i p (BioSequence t)) (BioSequence s) (BioSequence t)
{-# Inline pisSequence #-}
pisSequence = lens f t where
  v = view (locSequence.bioSequence)
  f (PIS p i s) = BioSequence $ maybe BS.empty v p `BS.append` v i `BS.append` maybe BS.empty v s
  t (PIS p i s) (BioSequence str) =
    let (pfx,ifxsfx) = over _1 BioSequence   $ BS.splitAt (maybe 0 (BS.length . v) p) str
        (ifx,sfx   ) = over both BioSequence $ BS.splitAt (BS.length $ v i) ifxsfx
    in  PIS (set (_Just . locSequence) pfx p) (set locSequence ifx i) (set (_Just . locSequence) sfx s)



-- | Given a @Location@ with a @BioSequence@, replace the sequence with its length.

locAsLength :: Location i FwdPosition (BioSequence w) -> Location i FwdPosition Int
{-# Inline locAsLength #-}
locAsLength = over locSequence (view (_BioSequence.to BS.length))



-- | Provides a range in a notation as used by blast, for example. This
-- isomorphism can translate back as well. @FwdLocation - 8 4 ^. blastRange1 ==
-- 9 6 MinusStrand@, since these ranges are 1-based and start and end included.

blastRange1 :: (Location i FwdPosition Int) -> (Int,Int,Strand)
{-# Inline blastRange1 #-}
blastRange1 = f -- iso f t
  where
    f loc =
      let s = loc^.locPosition.fwdStart.to toInt1
          l = loc^.locSequence
          pm = loc^.locPosition.fwdStrand
      in  case pm of PlusStrand -> (s,s+l,pm) ; MinusStrand -> (s+l,s,pm)
--    t (x,y,pm) =
--      let s = fromInt1 x
--          l = 1 + abs (x-y)
--      in  Location (FwdPosition pm s) l



-- | For each element, attach the prefix as well. The @Int@ indicates the maximal prefix length to
-- attach.
--
-- @1 2 3 4@ -> @01 12 23 34@
--
-- TODO are we sure this is correct for @MinusStrand@?

attachPrefixes
  :: ( Monad m, ModifyLocation p s )
  => Int
  -> SP.Stream (SP.Of (PIS i p s)) m r
  -> SP.Stream (SP.Of (PIS i p s)) m r
{-# Inlinable attachPrefixes #-}
attachPrefixes k = SP.map (\(Just w) -> w) . SP.drop 1 . SP.scan go Nothing id
  where
    go Nothing = Just
    go (Just p) = Just . set pisPrefix (Just . locTakeEnd k $ view pisInfix p)



-- | For each element, attach the suffix as well.
--
-- @1 2 3 4@ -> @12 23 34 40@

attachSuffixes
  :: ( Monad m, ModifyLocation p s )
  => Int
  -> SP.Stream (SP.Of (PIS i p s)) m r
  -> SP.Stream (SP.Of (PIS i p s)) m r
{-# Inlinable attachSuffixes #-}
attachSuffixes k = loop Nothing
  where
    loop Nothing = \case
      SI.Return r -> SI.Return r
      SI.Effect m -> SI.Effect $ fmap (loop Nothing) m
      SI.Step (a SP.:> rest) -> loop (Just a) rest
    loop (Just p) = \case
      SI.Return r -> SI.Step (p SP.:> SI.Return r)
      SI.Effect m -> SI.Effect $ fmap (loop (Just p)) m
      SI.Step (a SP.:> rest) ->
        let p' = p & set pisSuffix (Just . locTake k $ view pisInfix a)
        in  SI.Step (p' SP.:> loop (Just a) rest)