packages feed

bio 0.5.0.1 → 0.5.1

raw patch · 10 files changed

+182/−25 lines, 10 filesdep ~basenew-component:exe:flowclipnew-component:exe:orfnew-component:exe:rselect-pe

Dependency ranges changed: base

Files

Bio/Alignment/AlignData.hs view
@@ -12,7 +12,7 @@  -} -{-# OPTIONS -fglasgow-exts #-}+{-# LANGUAGE ParallelListComp #-}  module Bio.Alignment.AlignData (     -- * Data types for gap-based alignemnts
Bio/Sequence/FastQ.hs view
@@ -80,10 +80,13 @@ parse :: [B.ByteString] -> Maybe (Either String (Sequence Nuc), [B.ByteString]) parse (h1:sd:h2:sq:rest) =     case (B.uncons h1,B.uncons h2) of+      -- The fast path: four-line format       (Just ('@',h1name), Just ('+',h2name))-          | h1name == h2name -> Just (Right $ Seq h1name sd (Just (BB.map (subtract 33) sq)), rest)-          | otherwise        -> Just (Left $ "Bio.Sequence.FastQ: name mismatch:" ++ showStanza, rest)-      _                      -> Just (Left $ "Bio.Sequence.FastQ: illegal FastQ format:" ++ showStanza, rest)+          | h1name == h2name || B.null h2name+            -> Just (Right $ Seq h1name sd (Just (BB.map (subtract 33) sq)), rest)+          | otherwise+            -> Just (Left $ "Bio.Sequence.FastQ: name mismatch:" ++ showStanza, rest)+      _     -> Just (Left $ "Bio.Sequence.FastQ: illegal FastQ format:" ++ showStanza, rest)     where showStanza = unlines $ map B.unpack [ h1, sd, h2, sq ] parse [] = Nothing parse fs = let showStanza = unlines (map B.unpack fs)
Bio/Sequence/HashWord.lhs view
@@ -4,6 +4,8 @@    (Actually, it mostly packs nucleotide words from into integers, so it's a one-to-one relationship.)    It is intended to be used for indexing large sequence collections. +   Asking for a word size larger than will fit in the hash datatype will cause an error. Sorry.+ \begin{code} module Bio.Sequence.HashWord where @@ -40,6 +42,7 @@ contigous :: Integral k => Int -> HashF k contigous k' = let     k = fromIntegral k'+   klim = let v = 4^(k-1) in if v==0 then error ("HashWord: word size of k="++show k++" too large for this data type.") else v     c_key s i | k+i > B.length s = Nothing               | otherwise = let s' = B.take k $ B.drop i s @@ -55,7 +58,7 @@     c_keys'' cur p s | B.null s = []                     | isN (B.head s) = c_keys' (p+k+1) (B.tail s)-                    | otherwise = let new = (cur `mod` 4^(k-1))*4+val (B.head s)+                    | otherwise = let new = (cur `mod` klim)*4+val (B.head s)                                   in (new,p+1) : c_keys'' new (p+1) (B.tail s)                in HF { hash = c_key@@ -71,6 +74,7 @@ rcontig :: Integral k => Int -> HashF k rcontig k' = let    k = fromIntegral k'+   klim = let v = 4^(k-1) in if v==0 then error ("HashWord: word size of k="++show k++" too large for this data type.") else v     c_key s i | k+i > B.length s = Nothing               | B.find isN s' /= Nothing = Nothing@@ -90,8 +94,8 @@     c_keys'' (c1,c2) p s | B.null s = []                         | isN (B.head s) = c_keys' (p+k+1) (B.tail s)-                        | otherwise = let n1 = (c1 `mod` 4^(k-1))*4+val (B.head s)-                                          n2 = c2 `div` 4+4^(k-1)*(val . complement) (B.head s)+                        | otherwise = let n1 = (c1 `mod` klim)*4+val (B.head s)+                                          n2 = c2 `div` 4+klim*(val . complement) (B.head s)                                       in (min n1 n2,p+1) : c_keys'' (n1,n2) (p+1) (B.tail s)                in HF { hash = c_key@@ -121,6 +125,7 @@ rcpacked :: Integral k => Int -> HashF k rcpacked k' = let     k = fromIntegral k'+    klim = let v = 4^(k-1) in if v==0 then error ("HashWord: word size of k="++show k++" too large for this data type.") else v      c_key s i | k > B.length s' = Nothing               | B.any isN s'    = Nothing@@ -142,8 +147,8 @@     c_keys'' (c1,c2) i (s:ss)          | isN $ B.head s = c_keys' (i + (sum . map B.length) (s:take k' ss)) ss         | otherwise = let s1 = B.head s-                          n1 = (c1 `mod` 4^(k-1))*4+val s1-                          n2 = c2 `div` 4+4^(k-1)*(val . complement) s1+                          n1 = (c1 `mod` klim)*4+val s1+                          n2 = c2 `div` 4+klim*(val . complement) s1                       in (min n1 n2,i) : c_keys'' (n1,n2) (i+B.length s) ss               in HF { hash = c_key @@ -193,8 +198,8 @@ k2n k = fromStr . k2n' k  k2n' k i = if k==1 then [unval i]-           else let (q,r) = i `divMod` (4^(k-1)) in unval q : k2n' (k-1) r-+           else let (q,r) = i `divMod` klim in unval q : k2n' (k-1) r+             where klim = let v = 4^(k-1) in if v==0 then error ("HashWord: word size of k="++show k++" too large for this data type.") else v {- shiftRight, shiftLeft :: Integral k => Int -> Int -> k -> String -> k shiftRight k l pk str = (pk `mod` 4^l)*4^(k-l) + n2k 0 str
Bio/Sequence/Phd.hs view
@@ -30,11 +30,12 @@ mkPhd inp =    let (hd:fs) = filter (not . B.null) . B.lines $ inp       (comment,sd) = break (==B.pack "BEGIN_DNA") fs+      (sd', td) = break (==B.pack "END_DNA") sd       (magic,label) = B.splitAt 15 hd       more_magic = magic == B.pack "BEGIN_SEQUENCE "       fields = B.words . B.unlines                . filter (not . isSubstr (B.pack "_COMMENT")) $ comment-      sdata = filter ((==3).length) . map B.words $ sd+      sdata = filter ((==3).length) . map B.words $ sd'       err = error "failed to parse quality value"       qual = BB.fromChunks [BBB.pack . map (maybe err (fromIntegral . fst) . B.readInt . (!!1)) $ sdata]   in if more_magic then qual `seq` (Seq (compact $ B.unwords (label:fields))
Bio/Sequence/SFF.hs view
@@ -104,18 +104,16 @@  -- trimming the flowgram is necessary, but how to deal with the shift in flow -- sequence - i.e. what to do when trimming "splits" a flow into trimmed/untrimmed bases?--- | Trim a read to specific sequence position, inclusive bounds--- The current implementation has the unintended side effect of --- always trimming the flowgram down to a basecalled position.  --- Note that you can't (easily) write trimmed 'ReadBlock's to a file,--- since they need to have the same number of flows as given in --- the 'CommmonHeader'.++-- | Trim a read to specific sequence position, inclusive bounds. trimFromTo :: (Integral i) => i -> i -> ReadBlock -> ReadBlock trimFromTo x r rd = let   l = x-1   trim_seq = LB.drop (fromIntegral l) . LB.take (fromIntegral r)   trim_seq' = B.drop (fromIntegral l) . B.take (fromIntegral r)   trim_flw = B.drop ((2*) $ fromIntegral $ baseToFlowPos rd l) . B.take ((2*) $ fromIntegral $ baseToFlowPos rd r)+  new_flw  = trim_flw (flow_data rd)+  padding = B.replicate (B.length (flow_data rd) - B.length new_flw) 0   rh = read_header rd   [r',l'] = map fromIntegral [r,l]   rh' = rh { num_bases = fromIntegral (r'-l')@@ -123,7 +121,7 @@            , clip_qual_right = min (clip_qual_right rh-l') (r'-l'+1)            }   in rd { read_header = rh'-        , flow_data = trim_flw (flow_data rd)+        , flow_data = B.concat [new_flw, padding]         , flow_index = trim_seq' (flow_index rd)         , bases = trim_seq (bases rd)         , quality = trim_seq (quality rd)
Bio/Sequence/SeqData.hs view
@@ -108,9 +108,9 @@ seqToStr :: Sequence a -> Int -> Int -> [(Int, Int)] -> [Char] seqToStr (Seq header sd _) chunks len parts =     let strSeq     = toStr sd-        strHeader  = toStr header-        identifier = head (words strHeader)-        comment    = intercalate " " $ tail (words strHeader)+        strHeader  = words $ toStr header+        identifier = case strHeader of (i:_) -> i; _ -> ""+        comment    = intercalate " " $ case strHeader of (_:rs) -> rs; _ -> []          numLength  = (length . show . length) strSeq         width      = chunks * len + numLength + chunks
bio.cabal view
@@ -1,5 +1,5 @@ Name:                bio-Version:             0.5.0.1+Version:             0.5.1 License:             LGPL License-file:        LICENSE Cabal-Version:       >= 1.6@@ -39,7 +39,7 @@  Library   Build-Depends:     base>=4 && <5, binary >=0.4 && <0.5, tagsoup>=0.8, bytestring >= 0.9.1,-                     containers, array, parallel, parsec, mtl, directory+                     containers, array, parallel, parsec, mtl, directory, QuickCheck >= 2    Exposed-modules:   Bio.Sequence,                      Bio.Sequence.SeqData,@@ -68,7 +68,10 @@                      Bio.GFF3.SGD    Extensions:        CPP, ParallelListComp-  Ghc-Options:       -Wall -O2 -fexcess-precision -funbox-strict-fields -auto-all -fno-warn-unused-do-bind+  if impl(ghc >= 6.12)+    Ghc-Options:       -Wall -O2 -fexcess-precision -funbox-strict-fields -auto-all -fno-warn-unused-do-bind+  else+    Ghc-Options:       -Wall -O2 -fexcess-precision -funbox-strict-fields -auto-all  Source-Repository head   Type:           darcs@@ -126,6 +129,36 @@   Main-Is:        FRename.hs   Hs-Source-Dirs: examples .   Build-Depends:  base >= 3 && < 5, bytestring >= 0.9.1+  Ghc-Options:    -Wall+  if flag(examples)+     Buildable:   True+  else+     Buildable:   False++Executable     flowclip+  Main-Is:        FlowClip.hs+  Hs-Source-Dirs: examples .+  Build-Depends:  base >= 3 && < 5+  Ghc-Options:    -Wall+  if flag(examples)+     Buildable:   True+  else+     Buildable:   False++Executable     orf+  Main-Is:        Orf.hs+  Hs-Source-Dirs: examples .+  Build-Depends:  base >= 3 && < 5+  Ghc-Options:    -Wall+  if flag(examples)+     Buildable:   True+  else+     Buildable:   False++Executable     rselect-pe+  Main-Is:        RSelectPE.hs+  Hs-Source-Dirs: examples .+  Build-Depends:  base >= 3 && < 5, random   Ghc-Options:    -Wall   if flag(examples)      Buildable:   True
+ examples/FlowClip.hs view
@@ -0,0 +1,27 @@+{-| FlowClip reads in an SFF file using coordinates provided on the+    command line, and clips each read accordingly.  Negative coordinates+    indicate counting from the left, so to clip the last 20 bases off+    each read (useful for removing adapters), do 'flowclip 0 -20 input.sff'.+-}++module Main where++import Bio.Sequence.SFF+import System.Environment (getArgs)+import Data.Int++main :: IO ()+main = do+  [from,to,input] <- getArgs+  SFF h rs <- readSFF input+  writeSFF (input++".clip") (SFF h $ map (clip (read from) (read to)) rs)+  +clip :: Int32 -> Int32 -> ReadBlock -> ReadBlock+clip f t rb = let +  rh    = read_header rb+  n     = num_bases rh+  left  = if f < 0 then n+f else f+  right = if t < 0 then n+t else t+  in trimFromTo left right rb+  +
+ examples/Orf.hs view
@@ -0,0 +1,52 @@+-- find longest ORF in a Fasta file, inspired by a question on BioStar +-- (http://biostar.stackexchange.com/questions/5902/code-golf-finding-orf-and-corresponding-strand-in-a-dna-sequence/)++import Bio.Sequence+import System.IO+import Data.List (maximumBy, tails)++main :: IO ()+main = hReadFasta stdin >>= mapM_ (putStrLn . doit . castToNuc)++doit s = format s . maximumBy comp_orflength . all_orfs $ s++format :: Sequence Nuc -> (Frame,Length,Offset) -> String+format s (frm, len, off)+  | frm > 0 = sq++"+"++show frm++" "++show (frm+off)++" "++show (frm+off+len'-1)+  | otherwise = sq++show frm++" "++show (sl+frm-off-len'+2)++" "++show (sl-off+frm+1)+    where len' = fromIntegral len+          sl = fromIntegral (seqlength s)+          sq = toStr (seqheader s) ++ ":\t"++comp_orflength (_,l1,_) (_,l2,_) = compare l1 l2++type Frame = Offset -- -3..+3+type Length = Int++-- | Put it all togehter, and generate all ORFs from a sequence.+all_orfs :: Sequence Nuc -> [(Frame,Length,Offset)]+all_orfs = map (findlength . truncateStp) . concat . map transtails . translations++-- | Convert the translated sequence into the equivalent nucleotide length+findlength :: (Frame,Offset,[Amino]) -> (Frame,Length,Offset)+findlength (f,o,as) = (f,3*length as,o)++-- | Truncate the translation if and when a stop codon is reached.+truncateStp :: (Frame,Offset,[Amino]) -> (Frame,Offset,[Amino])+truncateStp (f,o,as) = (f,o,takeUntil (==STP) as)++-- | Just like takeWhile, but keeps the final element (ie. STP)+takeUntil p (x:xs) | p x = [x]+                   | otherwise = x : takeUntil p xs+takeUntil _ [] = []++-- | Analogous to the 'tails' function, this generates all suffixes+--   of a translation.+transtails :: (Frame,[Amino]) -> [(Frame,Offset,[Amino])]+transtails (f,as) = [(f,i,ts) | (i,ts) <- zip [0,3..] (tails as)]++-- | Generate all six-frame translations of the sequence, recording +--   the frame as well.  The order decides which one gets selected later        +--   on in the case of equally long ORFs.        +translations :: Sequence Nuc -> [(Frame,[Amino])]+translations s = [(negate i,translate (revcompl s) (i-1)) | i <- [3,2,1]] ++ [(i,translate s (i-1)) | i <- [3,2,1]]
+ examples/RSelectPE.hs view
@@ -0,0 +1,38 @@+{- | Randomly select paired end (illumina) reads -}++import System.Environment (getArgs)+import Data.List (isSuffixOf)+import System.Random+import System.IO+import Bio.Sequence++main = do+  (n:as) <- getArgs+  let fs = pairs sufchk as+      p  = read n :: Double+  rs <- map (<p) `fmap` randomRs (0,1.0) `fmap` newStdGen+  s1 <- concat `fmap` mapM readIllumina (map fst fs)+  s2 <- concat `fmap` mapM readIllumina (map snd fs)+  let out = [(s,t) | (s,t,True) <- zip3 s1 s2 rs]+  h1 <- openFile "out.1.txt" WriteMode +  h2 <- openFile "out.2.txt" WriteMode +  mapM_ (writer h1 h2) out+  hClose h1+  hClose h2+  +writer :: Handle -> Handle -> (Sequence Nuc,Sequence Nuc) -> IO ()+writer h1 h2 (s,t) | hdchk s t = do+  hWriteIllumina h1 [s]+  hWriteIllumina h2 [t]++pairs :: Show a => (a -> a -> Bool) -> [a] -> [(a,a)]  +pairs check (x:y:zs) +  | check x y = (x,y):pairs check zs+  | otherwise = error ("Suffix check failed for "++show x++" and "++show y)+pairs _ [] = []+pairs _ _x = error ("Uneven number of args, last was" ++ show _x)++sufchk x y = ".1.txt" `isSuffixOf` x && ".2.txt" `isSuffixOf` y++hdchk x y = f x == f y+  where f = takeWhile (/='/') . toStr . seqlabel