bio 0.3.3 → 0.3.3.1
raw patch · 7 files changed
+89/−55 lines, 7 filesdep ~tagsoupnew-uploader
Dependency ranges changed: tagsoup
Files
- Bio/Alignment/AAlign.hs +3/−3
- Bio/Alignment/ACE.hs +28/−33
- Bio/Alignment/AlignData.hs +44/−5
- Bio/Alignment/Multiple.hs +2/−2
- Bio/Alignment/QAlign.hs +6/−6
- Bio/Alignment/SAlign.hs +4/−4
- bio.cabal +2/−2
Bio/Alignment/AAlign.hs view
@@ -51,12 +51,12 @@ fp (x,ax) (s,e) = (x+s,e:ax) -- | Calculate global alignment (Needleman-Wunsch)-global_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,Alignment)+global_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,EditList) global_align mx g s1 s2 = revsnd . uncurry max' . last . last $ columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2 -- | Calculate local alignmnet (Smith-Waterman)-local_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,Alignment)+local_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,EditList) local_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') . concat $ columns (align_select 0 mx g) ((0,[]),(fst g,[])) s1 s2@@ -64,7 +64,7 @@ revsnd (s,a) = (s,reverse a) -- | Generic scoring and selection for global and local alignment-align_select :: (Num a, Ord a) => a -> SubstMx a -> (a,a) -> Selector ((a,Alignment),(a,Alignment))+align_select :: (Num a, Ord a) => a -> SubstMx a -> (a,a) -> Selector ((a,EditList),(a,EditList)) align_select minf mx (go,ge) cds = let (reps,ids) = partition (isRepl.snd) cds s = maximumBy (compare `on` fst)
Bio/Alignment/ACE.hs view
@@ -26,17 +26,15 @@ {-# LANGUAGE CPP #-} -module Bio.Alignment.ACE (readACE, writeACE, Assembly(..)- ,Gaps,extractGaps, insertGaps, ptest- )-where+module Bio.Alignment.ACE (readACE, writeACE, Assembly(..), ptest, reads ) where #include <interlude.h> import Interlude hiding (lines,words,readFile,unwords -- ByteString clashes- ,reads -- Assembly clash+ ,reads,pred -- Assembly clash ) import Bio.Sequence.SeqData+import Bio.Alignment.AlignData (Dir(..), Gaps, Alignment, extractGaps) import qualified Data.ByteString.Lazy.Char8 as B import Data.ByteString.Lazy.Char8 (ByteString,words,pack,unpack,readFile,unwords)@@ -47,11 +45,13 @@ import Control.Monad (liftM) -- ,when? import Data.Char (chr) -data Dir = Fwd | Rev deriving (Eq,Show)-data Assembly = Asm { contig :: (Sequence,Gaps), reads :: [(Offset,Dir,Sequence,Gaps)] }+data Assembly = Asm { contig :: (Sequence,Gaps), fragments :: Alignment } deriving Show -type Gaps = [Offset]+{-# DEPRECATED reads "Stupid name, replaced by 'fragments'." #-}+reads :: Assembly -> Alignment+reads = fragments -- deprecated, stupid, stupid name.+ type Str = ByteString -- | ACE header lines with parameters@@ -76,9 +76,9 @@ show Empty = "(blank)" show _ = "unknown ACE string" +uw :: [Str] -> String uw = unpack . unwords - -- | The Parsec parser type type AceParser a = GenParser (SourcePos,ACE) () a @@ -114,21 +114,10 @@ return (fst c,fst r) _ -> Nothing) <?> "AS <int> <int>" +blank :: AceParser () blank = parse1 (\t -> case t of Empty -> Just () _ -> Nothing) <?> "empty line" --- | Gaps are coded as '*'s, this function removes them, and returns--- the sequence along with the list of gap positions.-extractGaps :: Str -> (Str,Gaps)-extractGaps str = (B.filter (/='*') str,B.elemIndices '*' str)---- todo: faster to lift concat out of the inner loop?-insertGaps :: Char -> (Str,Gaps) -> Str-insertGaps c (str,gaps) = go str B.empty 0 gaps- where go str acc p (next:rest) = let (a,b) = B.splitAt (next-p) str- in go b (B.concat [acc,a,pack [c]]) (next+1) rest- go str acc _ [] = B.append acc str- -- | parse the contig and quality information (CO, BQ) ctg :: AceParser (Sequence,Gaps) ctg = do@@ -161,9 +150,9 @@ -- | Read a list of Ints in the Maybe monad readInts :: [ByteString] -> Maybe [Int] readInts [] = Just []-readInts (x:xs) = do (a,_) <- B.readInt x- as <- readInts xs- return $ (a:as)+readInts (x:xs) = do (i,_) <- B.readInt x+ is <- readInts xs+ return $ (i:is) bq :: AceParser () bq = parse1 (\t -> case t of BQ -> Just (); _ -> Nothing) <?> "BQ"@@ -179,24 +168,30 @@ -- | Parse a list of AFS, followed by actual read, and merge them -- afs :: Sequence -> AceParser [Sequence] -- plus some auxiliary info?--- todo: better error handling! af :: AceParser (Str,Dir,Offset)-af = parse1 (\t -> case t of AF a b c -> Just (a,f b,readInt' c)- _ -> Nothing) <?> "AF name compl pad_start"- where f b = case unpack b of "U" -> Fwd; "C" -> Rev- readInt' x = case B.readInt x of Just (a,_) -> fromIntegral a+af = parse1 (\t -> case t of AF a b c -> mkAF a b c+ _ -> Nothing) <?> "AF name (U|C) pad_start"+ where mkAF a b c = do+ b' <- case unpack b of "U" -> Just Fwd; "C" -> Just Rev; _ -> Nothing+ c' <- liftM (fromIntegral . fst) (B.readInt c)+ return (a,b',c') bs :: AceParser (Int,Int,Str)-bs = parse1 (\t -> case t of BS [x,y,n] -> Just (readInt' x, readInt' y,n)+bs = parse1 (\t -> case t of BS [x,y,n] -> do + x' <- readInt' x+ y' <- readInt' y+ return (x',y',n) _ -> Nothing) <?> "BS x y name"- where readInt' i = case B.readInt i of Just (a,_) -> fromIntegral a +readInt' :: Str -> Maybe Int+readInt' = liftM (fromIntegral . fst) . B.readInt+ rds :: (Sequence,Gaps) -> [(Str,Dir,Offset)] -> AceParser Assembly rds cg xs = do r <- many1 rseq -- todo: check the number and merge with the afs let f (_name,d,off) (s,gs) = (off,d,s,gs)- return $ Asm { contig = cg, reads = zipWith f xs r }+ return $ Asm { contig = cg, fragments = zipWith f xs r } rseq :: AceParser (Sequence,Gaps) rseq = do@@ -234,7 +229,7 @@ | otherwise = let (h:ws) = words l in case (unpack h,ws) of ("AS",[cs,rs]) -> AS cs rs- ("CO",[nm,bs,rs,segs,compl]) -> CO nm bs rs segs compl+ ("CO",[nm,bss,rs,segs,comp]) -> CO nm bss rs segs comp ("BQ",[]) -> BQ ("AF",[a,b,c]) -> AF a b c ("BS",_) -> BS ws
Bio/Alignment/AlignData.hs view
@@ -1,19 +1,58 @@--- | Data structures and helper functions for calculating alignments+{- | Data structures and helper functions for calculating alignments + There are two ways to view an alignment: either as a list of edits+ (i.e., insertions, deletions, or substitutions), or as a set of sequences+ with inserted gaps.++ The edit list approach is perhaps more restrictive model but doesn't generalize+ to multiple alignments.++ The gap approach is more general, and probably more commonly used by other software+ (see e.g. the ACE file format).++-}++{-# OPTIONS -fglasgow-exts #-}+ module Bio.Alignment.AlignData (- -- * Data types for alignments- Edit(..), Alignment, SubstMx, Selector, Chr+ -- * Data types for gap-based alignemnts+ Dir(..), Gaps, Alignment -- * Helper functions+ , extractGaps, insertGaps+ -- * Data types for edit-based alignments+ , Edit(..), EditList, SubstMx, Selector, Chr+ -- * Helper functions , columns, eval, isRepl, on , showalign, toStrings ) where import qualified Data.ByteString.Lazy as B+import qualified Data.ByteString.Lazy.Char8 as BC import Bio.Sequence.SeqData import Data.List (unfoldr) import Data.Word import Data.Char (chr) +-- ----------------------------------------++data Dir = Fwd | Rev deriving (Eq,Show)+type Gaps = [Offset]+type Alignment = [(Offset,Dir,Sequence,Gaps)]++-- | Gaps are coded as '*'s, this function removes them, and returns+-- the sequence along with the list of gap positions.+extractGaps :: SeqData -> (SeqData,Gaps)+extractGaps str = (BC.filter (/='*') str,BC.elemIndices '*' str)++-- todo: faster to lift concat out of the inner loop?+insertGaps :: Char -> (SeqData,Gaps) -> SeqData+insertGaps c (str',gaps) = go str' B.empty 0 gaps+ where go str acc p (next:rest) = let (a,b) = BC.splitAt (next-p) str+ in go b (BC.concat [acc,a,BC.pack [c]]) (next+1) rest+ go str acc _ [] = BC.append acc str++-- ----------------------------------------+ -- Q&D helper function showalign a = let (s1,s2) = toStrings a in s1++"\n"++s2 @@ -21,7 +60,7 @@ -- (for checking, filtering out the '-' characters should return -- the original sequences, provided '-' isn't part of the sequence -- alphabet)-toStrings :: Alignment -> (String,String)+toStrings :: EditList -> (String,String) toStrings [] = ("","") toStrings (x:xs) = let (a1',a2') = toStrings xs chr' = chr . fromIntegral@@ -37,7 +76,7 @@ data Edit = Ins Chr | Del Chr | Repl Chr Chr deriving (Show,Eq) -- | An alignment is a sequence of edits.-type Alignment = [Edit]+type EditList = [Edit] -- | True if the Edit is a Repl. isRepl :: Edit -> Bool
Bio/Alignment/Multiple.hs view
@@ -14,13 +14,13 @@ -- | Progressive multiple alignment. -- Calculate a tree from agglomerative clustering, then align -- at each branch going bottom up. Returns a list of columns (rows?).-progressive :: (Sequence -> Sequence -> (Double,Alignment)) -> [Sequence] -> [String]+progressive :: (Sequence -> Sequence -> (Double,EditList)) -> [Sequence] -> [String] progressive = undefined -- | Derive alignments indirectly, i.e. calculate A|C using alignments A|B and B|C. -- This is central for 'Coffee' evaluation of alignments, and T-Coffee construction -- of alignments.-indirect :: Alignment -> Alignment -> Alignment+indirect :: EditList -> EditList -> EditList indirect (Repl x1 x2:xs) (Repl y1 y2:ys) = Repl x1 y2 : indirect xs ys -- assert x2==y1 indirect xs@(Repl _ _:_) (Ins y1:ys) = Ins y1 : indirect xs ys indirect (Repl x1 _:xs) (Del y1:ys) = Del x1 : indirect xs ys
Bio/Alignment/QAlign.hs view
@@ -30,7 +30,7 @@ import Data.Array.Unboxed import Bio.Sequence.SeqData hiding ((!))-import Bio.Alignment.AlignData (Chr,Edit(..),SubstMx,Alignment,on,isRepl,showalign)+import Bio.Alignment.AlignData (Chr,Edit(..),SubstMx,EditList,on,isRepl,showalign) -- | The selector must take into account the quality of the sequences -- on Ins\/Del, the average of qualities surrounding the gap is (should be) used@@ -125,20 +125,20 @@ -- Alignments (rip from AAlign) -- maximum...-max' :: (Double,Alignment) -> (Double,Alignment) -> (Double,Alignment)+max' :: (Double,EditList) -> (Double,EditList) -> (Double,EditList) max' (x,ax) (y,yx) = if x>=y then (x,ax) else (y,yx) -- ... and addition for compound values-fp :: (Double,Alignment) -> (Double,Edit) -> (Double,Alignment)+fp :: (Double,EditList) -> (Double,Edit) -> (Double,EditList) fp (x,ax) (s,e) = (x+s,e:ax) -- | Calculate global alignment (Needleman-Wunsch)-global_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,Alignment)+global_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList) global_align mx g s1 s2 = revsnd . uncurry max' . last . last $ columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2 -- | Calculate local alignmnet (Smith-Waterman)-local_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double, Alignment)+local_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double, EditList) local_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') . concat $ columns (align_select 0 mx g) ((0,[]),(fst g,[])) s1 s2@@ -154,7 +154,7 @@ revsnd (s,a) = (s,reverse a) -- | Generic scoring and selection for global and local alignment-align_select :: Double -> QualMx Double -> (Double,Double) -> QSelector ((Double,Alignment),(Double,Alignment))+align_select :: Double -> QualMx Double -> (Double,Double) -> QSelector ((Double,EditList),(Double,EditList)) align_select minf mx (go,ge) cds = let (reps,ids) = partition (isRepl.snd') cds s = maximumBy (compare `on` fst)
Bio/Alignment/SAlign.hs view
@@ -46,17 +46,17 @@ -- the score in each cell. Unreachable cells can then be GC'ed. -- | Calculate alignments.-global_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> Alignment+global_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> EditList global_align mx g s1 s2 = reverse . snd . last . last $ columns (g_align mx g) (0,[]) s1 s2 -g_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,Alignment)+g_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,EditList) g_align mx g cds = maximumBy (compare `on` fst) [(s+eval mx g e,e:a) | ((s,a),e) <- cds] -local_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> Alignment +local_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> EditList local_align mx g s1 s2 = reverse . snd . maximumBy (compare `on` fst) . concat $ columns (l_align mx g) (0,[]) s1 s2 -l_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,Alignment)+l_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,EditList) l_align mx g cds = maximumBy (compare `on` fst) $ (0,[]):[(s+eval mx g e,e:a) | ((s,a),e) <- cds]
bio.cabal view
@@ -1,5 +1,5 @@ Name: bio-Version: 0.3.3+Version: 0.3.3.1 License: LGPL License-file: LICENSE Author: Ketil Malde@@ -24,7 +24,7 @@ Tested-With: GHC==6.8.2 Build-Type: Simple-Build-Depends: base>3, QuickCheck<2, binary, tagsoup>=0.5, bytestring,+Build-Depends: base>3, QuickCheck<2, binary, tagsoup>=0.4, bytestring, containers, array, interlude, parallel, parsec -- add fps for ghc 6.4.2; change imports in Bio/Sequence/TwoBit.hs if you want QC 2