packages feed

seqloc 0.6 → 0.6.1

raw patch · 4 files changed

+268/−8 lines, 4 files

Files

seqloc.cabal view
@@ -1,5 +1,5 @@ Name:                seqloc-Version:             0.6+Version:             0.6.1 Cabal-Version:       >= 1.14 Synopsis:            Handle sequence locations for bioinformatics Description:         Handle sequence locations for bioinformatics
src/Bio/SeqLoc/LocMap.hs view
@@ -1,28 +1,98 @@-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-module Bio.SeqLoc.LocMap+{-# LANGUAGE GeneralizedNewtypeDeriving, TypeSynonymInstances, FlexibleInstances #-}+{-| Efficient lookup based on potential location overlap++Collection of objects with locations, potentially on named sequences,+that can be queried to recover all objects whose location could+potentially overlap a query location.++|-}++module Bio.SeqLoc.LocMap (+  -- * Mapping objects to locations+  LocMap+  , emptyLM, insertLoc+  , queryLoc++  -- * Mapping objects to locations on named sequneces+  , SeqLocMap+  , emptySLM, insertSeqLoc+  , querySeqLoc++    -- * Mapping transcripts to their locations+  , transcriptSeqLocMap+    +  -- * Generalization of objects with locations+  , Locatable(..)+  , WithLocation(..)+  , locatableSeqLocMap+  , queryLocatable, queryLocCompatible+    +  )        where  import qualified Data.HashMap.Strict as HM import Data.List+-- import Data.Maybe  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 Bio.SeqLoc.Transcript  import qualified Bio.SeqLoc.ShiftedVector as ShV +-- | Mapping objects to sequence locations in the sense of `Loc.Location`. data LocMap a = LocMap { binSize :: !Pos.Offset,                          bins :: !(ShV.ShiftedVector [a]) } deriving (Show) +-- | Mapping objects to locations on named sequences  data SeqLocMap a = SeqLocMap { slmBinSize :: !Pos.Offset, locmaps :: !(HM.HashMap SeqLabel (LocMap a)) } deriving (Show) +-- | Object with a genomic location expressed by `SpliceSeqLoc`.+class Locatable o where+  locate :: o -> SpliceSeqLoc++instance Locatable ContigSeqLoc where+  locate (OnSeq ref loc) = OnSeq ref (SpLoc.singleton loc)++instance Locatable SpliceSeqLoc where+  locate = id++instance Locatable Transcript where+  locate = location++-- | Simple representation of a `Locatable` object as an arbitrary+-- object adjoined with a location.+data WithLocation a = WithLocation { withoutLocation :: !a, withLocate :: !SpliceSeqLoc } deriving (Eq)++instance Locatable (WithLocation a) where+  locate = withLocate++-- | Create an empty object / location map.+--+-- Specify a characteristic size for efficient queries, which depends+-- on the underlying genome. Smaller sizes yield fewer false+-- candidates but greater memory usage and potentially slower queries. emptyLM :: Pos.Offset -> LocMap a emptyLM bsz = LocMap { binSize = bsz, bins = ShV.empty } +-- | Insert a new object / location pair and reutrn the updated map insertLoc :: (Loc.Location l) => l -> a -> LocMap a -> LocMap a insertLoc l x lm0 = lm0 { bins = ShV.modifySome (bins lm0) (locBins l lm0) (x :)  } +-- | Retrieve a list of objects that could potentially overlap the+-- query location.+--+-- Some objects may not actually overlap the query location. Some+-- objects may appear more than once in the list. However, no object+-- whose location truly overlaps the query will be missing from the+-- list.+--+-- Overlap is defined by one or more nucleotides in common between the+-- bounds of the two locations, regardless of the relative strands of+-- the query and the object location. queryLoc :: (Loc.Location l) => l -> LocMap a -> [a] queryLoc l lm = concat [ (bins lm) ShV.!? b | b <- locBins l lm ] @@ -32,18 +102,77 @@                    binhigh = fromIntegral $ end `div` binSize lm                in [binlow..binhigh] +-- | Create an empty object / location map+--+-- Specify a characteristic size as per `emptyLM`. emptySLM :: Pos.Offset -> SeqLocMap a emptySLM bsz = SeqLocMap { slmBinSize = bsz, locmaps = HM.empty } +-- | Insert a new object / location pair and reutrn the updated map insertSeqLoc :: (Loc.Location l) => OnSeq l -> a -> SeqLocMap a -> SeqLocMap a insertSeqLoc sl x slm0 = let lm0 = HM.lookupDefault (emptyLM $ slmBinSize slm0) (onSeqLabel sl) (locmaps slm0)                              lm' = insertLoc (unOnSeq sl) x lm0                          in slm0 { locmaps = HM.insert (onSeqLabel sl) lm' (locmaps slm0) } +-- | Retrieve a list of objects that could potentially overlap the+-- query location.+--+-- Some objects may not actually overlap the query location. Some+-- objects may appear more than once in the list. However, no object+-- whose location truly overlaps the query will be missing from the+-- list.+--+-- Overlap is defined by one or more nucleotides in common between the+-- bounds of the two locations, regardless of the relative strands of+-- the query and the object location, as well as the same name for the+-- underlying reference sequence.+--+-- To retrieve objects and test for actual overlap see the `Locatable`+-- interface and `queryLocatable` or `queryLocCompatible`. querySeqLoc :: (Loc.Location l) => OnSeq l -> SeqLocMap a -> [a] querySeqLoc sl slm = maybe [] (queryLoc (unOnSeq sl)) $                      HM.lookup (onSeqLabel sl) (locmaps slm) +-- | Construct a mapping from transcripts to their genomic locations. transcriptSeqLocMap :: Pos.Offset -> [Transcript] -> SeqLocMap Transcript transcriptSeqLocMap bsz = foldl' insertTrx (emptySLM bsz)   where insertTrx slm0 t = insertSeqLoc (location t) t slm0++-- | Construct a mapping from a general `Locatable` object to its+-- genomic location+locatableSeqLocMap :: (Locatable l) => Pos.Offset -> [l] -> SeqLocMap l+locatableSeqLocMap bsz = foldl' insertLocable (emptySLM bsz)+  where insertLocable slm0 t = insertSeqLoc (locate t) t slm0++-- | Retrieve a list of `Locatable` objects whose overall, contiguous+-- genomic coordinates intersect at any position the genomic interval+-- spanned by the specified `Loc.Location`. This does not require that+-- the spliced structure of the query is a subset of the spliced+-- structure of the `Locatable` nor that the query location lie+-- entirely within the hit location (contrast with+-- `queryLocCompatible`).+--+-- When a strand argument is given, restrict matches to those lying on+-- the same strand as the query location, for `Just Plus`, or the+-- opposite strand, for `Just Minus`.+queryLocatable :: (Locatable o, Loc.Location l) => Maybe Strand -> OnSeq l -> SeqLocMap o -> [o]+queryLocatable mstr qyloc = filter realHit . querySeqLoc qyloc+  where realHit sb = let sbloc = locate sb+                         (qylb, qyub) = Loc.bounds . unOnSeq $ qyloc+                         (sblb, sbub) = Loc.bounds . unOnSeq $ sbloc+                         qystr = Loc.strand . unOnSeq $ qyloc+                         sbstr = Loc.strand . unOnSeq $ sbloc+                     in (qylb <= sbub) && (qyub >= sblb) && (qystr `strandCompat` sbstr)+        strandCompat = case mstr of+          Nothing -> \_ _ -> True+          Just Plus -> (==)+          Just Minus -> (/=)++-- | Retrieve a list of `Locatable` objects whose spliced structure+-- contains the query location specifically.+queryLocCompatible :: (Locatable o) => Maybe Strand -> SpliceSeqLoc -> SeqLocMap o -> [o]+queryLocCompatible mstr qyloc = filter compatHit . querySeqLoc qyloc+  where compatHit sb = maybe False strandCompat $ unOnSeq qyloc `SpLoc.locInto` (unOnSeq . locate $ sb)+        strandCompat = case mstr of+          Nothing -> const True+          Just s -> (== s) . Loc.strand
src/Bio/SeqLoc/SpliceLocation.hs view
@@ -11,8 +11,8 @@ module Bio.SeqLoc.SpliceLocation (    -- * Sequence locations   SpliceLoc-  , fromContigs-  , locOutof, locInto+  , fromContigs, singleton+  , locOutof, locInto, contigSpan   , mergeContigs, mergeAdjContigs   ) where  @@ -104,6 +104,11 @@ locInto :: (Location l) => SpliceLoc -> l -> Maybe SpliceLoc locInto sploc outer = mapM (flip clocInto outer) (toContigs sploc) >>=                       fromContigs . mergeContigs . concat . map toContigs++-- | Return the contiguous location spanned by the spliced location.+contigSpan :: SpliceLoc -> ContigLoc+contigSpan sl = let (start, end) = bounds sl+                in fromBoundsStrand start end (strand sl)  mergeContigs :: [ContigLoc] -> [ContigLoc] mergeContigs [] = []
test/TestMain.hs view
@@ -2,6 +2,7 @@ module Main     where +import Control.Applicative import Control.Monad import qualified Data.ByteString as BSW import qualified Data.ByteString.Char8 as BS@@ -43,6 +44,13 @@         , T "LocMap hit inside only"        property_LocMap_hitIn         , T "LocMap hit all"                property_LocMap_hitAll         , T "LocMap hit multi"              property_LocMap_hitMulti++        , T "Locatable hit Contig"           property_Locatable_hitContig+        , T "Locatable hit Spliced"          property_Locatable_hitSpliced+        , T "Locatable hit Contig Extended"  property_Locatable_locHitContig+        , T "Locatable hit Spliced Extended" property_Locatable_locHitSpliced+        , T "Locatable compat Contig Ext'd"  property_Locatable_locCompatContig+        , T "Locatable compat Spliced Ext'd" property_Locatable_locCompatSpliced                    , T "Pos revCompl"                  test_Pos_revCompl         , T "Pos atPos"                     property_Pos_atPos@@ -445,7 +453,7 @@       ploc = Loc.fromPosLen pos 1       lm = LM.insertLoc contig contig (LM.emptyLM binsz)       ishit = not $ null (LM.queryLoc ploc lm)-  in collect ishit $ ishit || not isin+  in ishit || not isin  property_LocMap_hitMulti :: Loc.ContigLoc -> Loc.ContigLoc -> Loc.ContigLoc -> Pos.Pos -> Property property_LocMap_hitMulti ca cb cc pos =@@ -460,10 +468,128 @@       hita = ca `elem` LM.queryLoc ploc lm       hitb = cb `elem` LM.queryLoc ploc lm       hitc = cc `elem` LM.queryLoc ploc lm-  in collect (hita, hitb, hitc) $-     and [ hita || not isina,+  in and [ hita || not isina,            hitb || not isinb,            hitc || not isinc ]++property_Locatable_hitContig :: [Loc.ContigLoc] -> Pos.Pos -> Property+property_Locatable_hitContig cs pos =+  forAll genPositiveOffset $ \binsz ->+  forAll genName $ \chr ->+  let isins = map (isJust . Loc.posInto pos) cs+      ploc = OnSeq chr $ Loc.fromPosLen pos 1+      seqlocs = map ((OnSeq chr) . SpLoc.singleton) cs+      locables = [ LM.WithLocation sl sl | sl <- seqlocs ]+      slm = LM.locatableSeqLocMap binsz locables+      binhits = map (\l -> l `elem` LM.querySeqLoc ploc slm) locables+      hits = map (\l -> l `elem` LM.queryLocatable Nothing ploc slm) locables+  in and $ concat [ zipWith (\isin binhit -> binhit || not isin) isins binhits+                  , zipWith (\isin hit -> hit == isin) isins hits+                  ]++property_Locatable_hitSpliced :: [SpLoc.SpliceLoc] -> Pos.Pos -> Property+property_Locatable_hitSpliced cs pos =+  forAll genPositiveOffset $ \binsz ->+  forAll genName $ \chr ->+  let isins = map (isJust . Loc.posInto pos) cs+      ploc = OnSeq chr . SpLoc.singleton $ Loc.fromPosLen pos 1+      seqlocs = map (OnSeq chr) cs+      locables = [ LM.WithLocation sl sl | sl <- seqlocs ]+      slm = LM.locatableSeqLocMap binsz locables+      binhits = map (\l -> l `elem` LM.querySeqLoc ploc slm) locables+      hits = map (\l -> l `elem` LM.queryLocatable Nothing ploc slm) locables+      compats = map (\l -> l `elem` LM.queryLocCompatible Nothing ploc slm) locables+  in and $ concat [ zipWith (\isin binhit -> binhit || not isin) isins binhits+                  , zipWith (\isin hit -> hit || not isin) isins hits+                  , zipWith (\isin compat -> isin == compat) isins compats+                  ]++property_Locatable_locHitContig :: [Loc.ContigLoc] -> Loc.ContigLoc -> Property+property_Locatable_locHitContig cs qy =+  forAll genPositiveOffset $ \binsz ->+  forAll genName $ \chr ->+  let isInFwds = map (Loc.overlaps qy) cs+      isInRevs = map (Loc.overlaps (revCompl qy)) cs+      isIns = zipWith (||) isInFwds isInRevs+      qyseqloc = OnSeq chr qy+      seqlocs = map ((OnSeq chr) . SpLoc.singleton) cs+      locables = [ LM.WithLocation sl sl | sl <- seqlocs ]+      slm = LM.locatableSeqLocMap binsz locables+      binhits = map (\l -> l `elem` LM.querySeqLoc qyseqloc slm) locables+      fwdHits = map (\l -> l `elem` LM.queryLocatable (Just Plus) qyseqloc slm) locables+      revHits = map (\l -> l `elem` LM.queryLocatable (Just Minus) qyseqloc slm) locables+      allHits = zipWith (||) fwdHits revHits+  in and $ concat [ zipWith (\isin binhit -> binhit || not isin) isIns binhits+                  , zipWith (\isin hit -> hit == isin) isInFwds fwdHits+                  , zipWith (\isin hit -> hit == isin) isInRevs revHits+                  , zipWith (\isin hit -> hit == isin) isIns allHits+                  ]++property_Locatable_locHitSpliced :: [SpLoc.SpliceLoc] -> Loc.ContigLoc -> Property+property_Locatable_locHitSpliced ls qy =+  forAll genPositiveOffset $ \binsz ->+  forAll genName $ \chr ->+  let isInFwds = map (Loc.overlaps qy . SpLoc.contigSpan) ls+      isInRevs = map (Loc.overlaps (revCompl qy) . SpLoc.contigSpan) ls+      isIns = zipWith (||) isInFwds isInRevs+      qyseqloc = OnSeq chr qy+      seqlocs = map (OnSeq chr) ls+      locables = [ LM.WithLocation sl sl | sl <- seqlocs ]+      slm = LM.locatableSeqLocMap binsz locables+      binhits = map (\l -> l `elem` LM.querySeqLoc qyseqloc slm) locables+      fwdHits = map (\l -> l `elem` LM.queryLocatable (Just Plus) qyseqloc slm) locables+      revHits = map (\l -> l `elem` LM.queryLocatable (Just Minus) qyseqloc slm) locables+      allHits = zipWith (||) fwdHits revHits+  in and $ concat [ zipWith (\isin binhit -> binhit || not isin) isIns binhits+                  , zipWith (\isin hit -> hit == isin) isInFwds fwdHits+                  , zipWith (\isin hit -> hit == isin) isInRevs revHits+                  , zipWith (\isin hit -> hit == isin) isIns allHits+                  ]++property_Locatable_locCompatContig :: [Loc.ContigLoc] -> Loc.ContigLoc -> Property+property_Locatable_locCompatContig cs qy =+  forAll genPositiveOffset $ \binsz ->+  forAll genName $ \chr ->+  let isInFwds = map ((== (Just Plus)) . liftM Loc.strand . Loc.clocInto qy) cs+      isInRevs = map ((== (Just Minus)) . liftM Loc.strand . Loc.clocInto qy) cs+      isIns = zipWith (||) isInFwds isInRevs+      qyseqloc = OnSeq chr (SpLoc.singleton qy)+      seqlocs = map ((OnSeq chr) . SpLoc.singleton) cs+      locables = [ LM.WithLocation sl sl | sl <- seqlocs ]+      slm = LM.locatableSeqLocMap binsz locables+      binhits = map (\l -> l `elem` LM.querySeqLoc qyseqloc slm) locables+      fwdHits = map (\l -> l `elem` LM.queryLocCompatible (Just Plus) qyseqloc slm) locables+      revHits = map (\l -> l `elem` LM.queryLocCompatible (Just Minus) qyseqloc slm) locables+      allHits = zipWith (||) fwdHits revHits+  in and $ concat [ zipWith (\isin binhit -> binhit || not isin) isIns binhits+                  , zipWith (\isin hit -> hit == isin) isInFwds fwdHits+                  , zipWith (\isin hit -> hit == isin) isInRevs revHits+                  , zipWith (\isin hit -> hit == isin) isIns allHits+                  ]++property_Locatable_locCompatSpliced :: [SpLoc.SpliceLoc] -> Loc.ContigLoc -> Property+property_Locatable_locCompatSpliced ls qy =+  forAll genPositiveOffset $ \binsz ->+  forAll genName $ \chr ->+  let isInFwds = map ((== (Just Plus)) . liftM Loc.strand . Loc.clocInto qy) ls+      isInRevs = map ((== (Just Minus)) . liftM Loc.strand . Loc.clocInto qy) ls+      isIns = zipWith (||) isInFwds isInRevs+      qyseqloc = OnSeq chr (SpLoc.singleton qy)+      seqlocs = map (OnSeq chr) ls+      locables = [ LM.WithLocation sl sl | sl <- seqlocs ]+      slm = LM.locatableSeqLocMap binsz locables+      binhits = map (\l -> l `elem` LM.querySeqLoc qyseqloc slm) locables+      fwdHits = map (\l -> l `elem` LM.queryLocCompatible (Just Plus) qyseqloc slm) locables+      revHits = map (\l -> l `elem` LM.queryLocCompatible (Just Minus) qyseqloc slm) locables+      allHits = zipWith (||) fwdHits revHits+  in and $ concat [ zipWith (\isin binhit -> binhit || not isin) isIns binhits+                  , zipWith (\isin hit -> hit == isin) isInFwds fwdHits+                  , zipWith (\isin hit -> hit == isin) isInRevs revHits+                  , zipWith (\isin hit -> hit == isin) isIns allHits+                  ]++instance (Arbitrary a) => Arbitrary (OnSeq a) where+  arbitrary = OnSeq <$> arbitrary <*> arbitrary  -- Utilities