trimpolya (empty) → 0

raw patch · 3 files changed

+68/−0 lines, 3 filesdep +basedep +biodep +bytestringsetup-changed

Dependencies added: base, bio, bytestring, simpleargs

Files

• Setup.hs +5/−0
• TrimPolyA.hs +39/−0
• trimpolya.cabal +24/−0

@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell+ +import Distribution.Simple+main = defaultMain+

@@ -0,0 +1,39 @@+module Main where++import qualified Data.ByteString.Lazy as BB+import qualified Data.ByteString.Lazy.Char8 as B+import Data.Char (toUpper)++import Bio.Sequence+import System.SimpleArgs++main = do+  (fs,qs) <- getArgs+  trimPolyA fs qs++trimPolyA fs qs = writeFastaQual (fs++".trim") (qs++".trim") =<< map doTrim `fmap` map castToNuc `fmap` readFastaQual fs qs++doTrim s = case findPolyA s of Just (x,y) -> slice (x,y) s+                               Nothing -> s++slice :: (Int,Int) -> Sequence Nuc -> Sequence Nuc+slice (x,y) (Seq h d (Just q)) = (Seq (B.concat [h,B.pack (" trim: "++show (y+1)++" poly-A: "++show (y-x+1))]) +                                      (B.take (fromIntegral y) d) (Just $ B.take (fromIntegral y) q))++-- Ripped from Dephd+findPolyA :: Sequence Nuc -> Maybe (Int,Int)+findPolyA (Seq _ d mq) = +      let qd = zip (B.unpack d) (maybe (repeat 15) BB.unpack mq)+          scores = map (\(c,q) -> if toUpper c=='A' then match q else mismatch q) qd+          match x' = let x = fromIntegral x' in log (4*(1-1/10**(x/10)))+          mismatch x' = let x = fromIntegral x' in log 4 - log 3 - x/10*log 10+          cumulative = scanl (\a b -> let r = a + b in max 0 r) 0+          (zi,mi,maxscore) = findmax $ cumulative scores+      in if maxscore > 12 then Just (zi+1,mi) else Nothing  -- arbitrary constant alert!++findmax :: [Double] -> (Int,Int,Double)+findmax = go 0 (0,0,0) . zip [0..]+    where go _ cm [] = cm+          go _ cm ((i,0):rest) = go i cm rest+          go last_z (cmz,cmi,cmx) ((i,x):rest) = if x > cmx then go last_z (last_z,i,x) rest +                                                 else go last_z (cmz,cmi,cmx) rest

@@ -0,0 +1,24 @@+Name:           trimpolya+Version:        0+License:        GPL++Author:         Ketil Malde+Maintainer:     Ketil Malde <ketil@malde.org>++Category:       Bioinformatics+Synopsis:       Search for, annotate and trim poly-A tail+Description:    trimpolya - a very simple tool for identifying and annotating poly-A tails+                and trimming following sequence. +                .+                Usage: trimpolya input.fasta input.qual+		.+                Detailed description of the method is at <http://blog.malde.org/index.php/2009/12/14/searching-for-polya-tails/>,+                darcs repository is at <http://malde.org/~ketil/biohaskell/trimpolya>.++Build-Depends:  base>=3 && <4, bio > 0.4, simpleargs, bytestring+Build-Type:     Simple+Tested-with:    GHC==6.8.3++Executable:     trimpolya+Main-Is:        TrimPolyA.hs+Ghc-Options:    -Wall