diff --git a/Bio/Alignment/AlignData.hs b/Bio/Alignment/AlignData.hs --- a/Bio/Alignment/AlignData.hs +++ b/Bio/Alignment/AlignData.hs @@ -12,7 +12,7 @@ -} -{-# OPTIONS -fglasgow-exts #-} +{-# LANGUAGE ParallelListComp #-} module Bio.Alignment.AlignData ( -- * Data types for gap-based alignemnts diff --git a/Bio/Sequence/FastQ.hs b/Bio/Sequence/FastQ.hs --- a/Bio/Sequence/FastQ.hs +++ b/Bio/Sequence/FastQ.hs @@ -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) diff --git a/Bio/Sequence/HashWord.lhs b/Bio/Sequence/HashWord.lhs --- a/Bio/Sequence/HashWord.lhs +++ b/Bio/Sequence/HashWord.lhs @@ -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 diff --git a/Bio/Sequence/Phd.hs b/Bio/Sequence/Phd.hs --- a/Bio/Sequence/Phd.hs +++ b/Bio/Sequence/Phd.hs @@ -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)) diff --git a/Bio/Sequence/SFF.hs b/Bio/Sequence/SFF.hs --- a/Bio/Sequence/SFF.hs +++ b/Bio/Sequence/SFF.hs @@ -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) diff --git a/Bio/Sequence/SeqData.hs b/Bio/Sequence/SeqData.hs --- a/Bio/Sequence/SeqData.hs +++ b/Bio/Sequence/SeqData.hs @@ -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 diff --git a/bio.cabal b/bio.cabal --- a/bio.cabal +++ b/bio.cabal @@ -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 diff --git a/examples/FlowClip.hs b/examples/FlowClip.hs new file mode 100644 --- /dev/null +++ b/examples/FlowClip.hs @@ -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 + + diff --git a/examples/Orf.hs b/examples/Orf.hs new file mode 100644 --- /dev/null +++ b/examples/Orf.hs @@ -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]] diff --git a/examples/RSelectPE.hs b/examples/RSelectPE.hs new file mode 100644 --- /dev/null +++ b/examples/RSelectPE.hs @@ -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 ( 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