diff --git a/seqloc.cabal b/seqloc.cabal --- a/seqloc.cabal +++ b/seqloc.cabal @@ -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 diff --git a/src/Bio/SeqLoc/LocMap.hs b/src/Bio/SeqLoc/LocMap.hs --- a/src/Bio/SeqLoc/LocMap.hs +++ b/src/Bio/SeqLoc/LocMap.hs @@ -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 diff --git a/src/Bio/SeqLoc/SpliceLocation.hs b/src/Bio/SeqLoc/SpliceLocation.hs --- a/src/Bio/SeqLoc/SpliceLocation.hs +++ b/src/Bio/SeqLoc/SpliceLocation.hs @@ -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 [] = [] diff --git a/test/TestMain.hs b/test/TestMain.hs --- a/test/TestMain.hs +++ b/test/TestMain.hs @@ -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