packages feed

dephd 0.1.3 → 0.1.4

raw patch · 4 files changed

+240/−34 lines, 4 filesdep +cmdargsdep ~basedep ~bionew-component:exe:fakequal

Dependencies added: cmdargs

Dependency ranges changed: base, bio

Files

README view
@@ -40,13 +40,21 @@ use -X on multiple files, hit q to terminate one window and go to the next. -The remaining three options (-R, -F, and -Q) output different aspects+Three other options (-R, -F, and -Q) output different aspects of the sequence information to files (specify '-' for printing to standard output instead - obviously this will be messy if you do it for more than one option!).  -F and -Q is for generating the standard Fasta and Quality files, while -R produces a file with one line per sequence containing various quality measures, including a verdict ranging from Excellent, through Good and Poor, to Junk.++Finally, -E can be used to generate a file suitable for submission to+dbEST.  Usually, you need to provide a library table (-l option), that is,+a text file with whitespace-seprated columns describing each library.+The first line of the table contains columns labels, and these should+correspond to fields in the dbEST library record format. One column+should be labelled "Pattern" and contain a regular expression matching+sequence names from this library.  Filtering can be specified with the -t option, which interprets trimming information from Phred or Lucy, and chops off the offending
dephd.cabal view
@@ -1,5 +1,5 @@ Name:           dephd-Version:        0.1.3+Version:        0.1.4 License:        GPL License-File:   LICENSE @@ -17,17 +17,28 @@ 		. 		Can trim according to Lucy or Phred parameters, can mask by quality, can plot 		graphs (via gnuplot) of sequence quality to a window, or to JPG/EPS files.  Can-		categorize sequences according to overall quality.+		categorize sequences according to overall quality.  Also constructs files suitable for+		submission to dbEST.  More information at <http://blog.malde.org/index.php/2010/09/07/submitting-ests-upstream/>. 		.+		Also provides 'fakequal', a utility to generate bogus quality values,+		which are sometimes needed by less flexible tools.+		.                 The Darcs repository is at <http://malde.org/~ketil/biohaskell/dephd>.  HomePage:	http://malde.org/~ketil/biohaskell/dephd-Build-Depends:  base>=3 && <4, bio >= 0.4, regex-compat, bytestring, process, directory+Build-Depends:  base>=3 && <5, bio > 0.4, regex-compat, bytestring, process, directory, cmdargs Build-Type:     Simple-Tested-with:    GHC==6.8.3+Cabal-Version:	>= 1.2.3+Tested-with:    GHC==6.12.1  Data-files:     README, TODO Executable:     dephd Main-Is:        Dephd.hs Hs-Source-Dirs: src+Extensions: 	BangPatterns Ghc-Options:    -Wall++Executable: 	fakequal+Main-Is:	FakeQual.hs+Hs-Source-Dirs: src+Extensions:	DeriveDataTypeable
src/Dephd.hs view
@@ -3,16 +3,18 @@    Generates seq/qual output, quality plots, or rankings. -} +{-# LANGUAGE BangPatterns #-}+ module Main where  import Control.Concurrent import Control.Monad import Data.Char-import Data.List (tails,groupBy,sortBy,isPrefixOf)+import Data.List (groupBy,sortBy,isPrefixOf,unfoldr) import Data.Maybe import System.Console.GetOpt import System.Directory-import System.Environment (getArgs)+import System.Environment (getArgs,getEnvironment) import System.Exit import System.IO import System.IO.Unsafe@@ -20,8 +22,10 @@ import Text.Printf import Text.Regex import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.ByteString.Lazy as BB  import Bio.Sequence+import Bio.Sequence.SeqData (hasqual) import Bio.Util (countIO)  -- ------------------------------------------------------------@@ -34,11 +38,13 @@                 , zerofilter :: (FilePath,Sequence Nuc) -> Bool                 , inputs  :: [String] -> IO [(FilePath,Sequence Nuc)] -- ^ Turn args into sequences                 , verbose :: Bool  -- ^ Verbose output (progress reporting) and sequence trimming+                , libtable :: Maybe FilePath                 }  defaultopts :: MyOpts defaultopts = O { actions = [], outputs = [], filters = id                 , zerofilter = const True+                , libtable = Nothing                 , inputs = readFiles, verbose = False }  getOptions :: IO (MyOpts, [String], [String])@@ -58,14 +64,23 @@     , Option ['R'] ["output-ranks"]  (ReqArg setrank "file") "Set ranked output file"     , Option ['F'] ["output-fasta"] (ReqArg setfasta "file") "Set fasta output file"     , Option ['Q'] ["output-qual"]  (ReqArg setqual "file")  "Set quality output file"+    , Option ['E'] ["output-dbEST"] +        (ReqArg (\arg opt -> do h <- openFile arg WriteMode +                                lt <- getLibTable (libtable opt) h+                                (c,p) <- getContPub h+                                w <- mkWriteEST lt (c,p)+                                return opt { actions = (w h.snd):actions opt, outputs = h:outputs opt}) "file")   "Output file suitable for dbEST submission"     , Option ['J'] ["output-plot"]  (NoArg  (setplot J))   "Plot sequence quality, JPEG"     , Option ['P'] ["output-plot"]  (NoArg  (setplot P))   "Plot sequence quality, EPS"     , Option ['X'] ["output-xplot"] (NoArg  (setplot X))   "Display quality plots" +    , Option ['l'] ["libtable"] (ReqArg (\arg opt -> return opt { libtable = Just arg }) "FILE")  "Specify a library table for ESTs"+     -- Filter options     , Option ['t'] ["filter-trim"] (NoArg filterTrim) "Trim output sequences.\nSpecify *before* -q if you want to trim based on quality!"-    , Option ['q'] ["filter-qual"] (NoArg filterQual) "Mask by quality"+    , Option ['q'] ["filter-qual"] (ReqArg filterQual "num") "Mask by quality to given threshold"     -- Input options+             , Option [] ["input-dirs"]     (NoArg  (\opt -> return opt { inputs = readDirs }))                                                             "Read directories containing PHD files"     , Option [] ["input-list"]     (NoArg inputList) "Read the files listed in an index file"@@ -84,7 +99,8 @@                                         return $ map (\s -> (toStr $ seqlabel s, castToNuc s)) ss }       addFilter act opt     = let f = filters opt in return $ opt { filters = f . act }       filterTrim = addFilter trimSeq-      filterQual = addFilter (\(p,s) -> (p,qualAdjust s))+      filterQual t = case reads t of [(t',_)] -> addFilter (\(p,s) -> (p,qualAdjust t' s))+                                     _ -> error ("Couldn't parse numeric argument to -q: '"++t++"'")       zeroFilter = (/=0). seqlength . snd       addAction act arg opt = let as = actions opt                                   hs = outputs opt@@ -101,7 +117,6 @@       setrank  = addAction (\h -> hPutStrLn h . unwords {- columns -} . qualchk . snd)       setfasta = addAction (\h -> hWriteFasta h . return . snd)       setqual  = addAction (\h -> hWriteQual h . return . snd)-       setplot c opt = do              g <- hasgp              let how = case c of J -> \f -> "set terminal jpeg\nset output \""@@ -112,6 +127,112 @@              if g then let as = actions opt in return $ opt { actions = plot how : as }                   else error ("You requested quality plots, but I can't find 'gnuplot' in your search path.\n") ++type LibTable = [[(String,String)]]  -- | column and value++readLibTable :: FilePath -> IO LibTable+readLibTable f = do+  let decomment = filter (\x -> (not . null) x && head x /= '#')+  (ch:ls) <- (map words . decomment . lines) `fmap` readFile f+  return $ map (zip ch) ls++-- | Classify a read according to the LibTable+classify :: [(Regex,String)] -> String -> String+classify ps str = case concatMap (class1 str) ps of+                    [] -> error ("no match for "++str++" in library table")+                    [x] -> x+                    s@(_:_) -> error ("multiple matches for "++str++": "++show s)+    where class1 st (r,s) = maybe [] (const [s]) (matchRegex r st)+++-- Write ESTs in (partial) GenBank/dbEST format+-- (todo: speed up by going to ByteStrings)++getContPub :: Handle -> IO (String,String)+getContPub h = do  +  let fs = ["contact.txt","publication.txt"]+  ftest <- mapM doesFileExist fs+  when (not $ and ftest) $ do+       mapM_ mkContPub $ map snd $ filter (not . fst) $ zip ftest fs+       error ("Please fill in the required information")+  let field x f = do ls <- (filter ((x++":") `isPrefixOf`) . lines) `fmap` readFile f +                     case ls of [l] -> return (dropWhile isSpace $ drop (length x+1) l)+                                _ -> error ("File: "++f++" is present, but does not contain field "++show x)+  readFile "publication.txt" >>= hPutStr h+  hPutStrLn h "||"+  readFile "contact.txt" >>= hPutStr h  +  hPutStrLn h "||"  +  +  cont_name <- field "NAME" "contact.txt"+  citation  <- field "TITLE" "publication.txt"+  return (cont_name,citation)  ++-- | Generate skeleton files if they don't exist.+-- contact.txt must contain a NAME field, "citation.txt" must contain TITLE.+mkContPub :: String -> IO ()+mkContPub f = do+  writeFile f $ maybe "" id $ lookup f skeletons+  putStrLn ("There was no file '"++f+    ++"', I generated one for you, please make sure it contains the appropriate data.")+  where skeletons = [("contact.txt","TYPE:\tCont\nNAME:\tobligatory!\nTEL:\nEMAIL:\nINST:\nADDR:\n")+                    ,("publication.txt","TYPE:\tPub\nTITLE:\tobligatory!\nAUTHORS:\nJOURNAL:\nVOLUME:\nISSUE:\nYEAR:\nSTATUS:\n")]++getLibTable :: Maybe FilePath -> Handle -> IO LibTable+getLibTable Nothing _ = error "Specify a libtable before -E"+getLibTable (Just lt) h = do+  libtab <- readLibTable lt+  let split str = let (k,v) = span (/=':') str+                  in if null k || null v then error ("incorrect format in $DBLIB: "++str)+                     else (k,dropWhile isSpace $ drop 1 v)+  extraLibs <- map split `fmap` maybe [] lines `fmap` lookup "DBLIB" `fmap` getEnvironment+  let trsp '_' = ' ' -- underscores represent spaces in libtable+      trsp x   = x+      extract fs n = case filter ((==(map toLower n)) . map toLower . fst) (fs++extraLibs) of+        [] -> if n `elem` ["ORGANISM","NAME"]+              then error (n++" is a required field - specify in "++lt++" or $DBLIB")+              else ""+        [(_,x)] -> if x/="?" then n++":\t"++map trsp x++"\n" else ""+        _ -> error ("Multiple definitions for field \""++n++"\"")+      genlib fs =  "TYPE:\tLib\n"+                   ++ concatMap (extract fs) ["NAME", "TISSUE", "STAGE", "ORGANISM", "STRAIN"]+                   ++ "||\n"+  hPutStr h $ concatMap genlib libtab+  return libtab+               +mkWriteEST :: LibTable -> (String,String) -> IO (Handle -> Sequence Nuc -> IO ())+mkWriteEST libtab  (cont_name,citation) = do+  let matchtable = let ns = concatMap (map snd . filter ((=="Name") . fst)) $ libtab+                       rs = concatMap (map snd . filter ((=="Pattern") . fst)) $ libtab+                   in [(mkRegex r,n) | (r,n) <- zip rs ns]++      lookuplibrary = classify matchtable++  extraFields <- maybe [] lines `fmap` lookup "DBEST" `fmap` getEnvironment++  let hwe h s = let hiqual = if not (hasqual s) then []+                               else let as = sliding_avg 20 (seqqual s)+                                        (trim_left,trim_right) = (length $ takeWhile (<20) as, fromIntegral (seqlength s) - length (takeWhile (<20) $ reverse as))+                               in ["HIQUAL_START:\t" ++ show (1+min trim_left trim_right), "HIQUAL_STOP:\t"  ++ show trim_right] -- if trim_right is less, there is no high qual+                    polya = [case findPolyA s of Just _ -> "POLYA:\tY"; Nothing -> "POLYA:\tN"]     -- POLYA: "Y" or "N"+                    clone   = []   -- ditto?+                    put_id  = []   -- from annotations.csv+                    est_name = toStr (seqlabel s) +                in case seqlength s of+                  0 -> return ()+                  _ -> hPutStr h $ unlines $+                         [ "TYPE:\tEST"+                         , "STATUS:\tNew"+                         , "CONT_NAME:\t"++cont_name+                         , "CITATION:\t"++citation+                         , "LIBRARY:\t"++lookuplibrary est_name+                         , "EST#:\t"++ est_name+                         ] ++ clone ++ put_id ++ hiqual ++ polya ++ extraFields +++                         [ "COMMENT:\tgenerated by dephd"+                         , "SEQUENCE:\t"++toStr (seqdata s)+                         , "||"+                         ]+  return hwe+ trimSeq :: (FilePath,Sequence Nuc) -> (FilePath,Sequence Nuc) trimSeq (i,s@(Seq _ d mq)) =      case trims s of @@ -125,6 +246,26 @@ -- todo: add trimming info in header?   -- (Currently, nothing protects agains re-trimming with the same parameters...) +-- PolyA finding.  Given error probability e, the probability that the base really is A is +-- match: 1-e/0.25   mismatch e/3/0.25.  Qual Q => e = 10^(-Q/10) => 1-10^(-Q/10)/0.25+-- ie. match: log 4 - Q/10*log 10 - log 3+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+ hasgp :: IO Bool hasgp = return . isJust =<< findExecutable "gnuplot" @@ -167,20 +308,26 @@ myGetDirectoryContents d = return . map ((d++"/")++) =<< getDirectoryContents d  -- | Adjust sequence content according to quality.---   Upper case is >20 and sliding avg >25-qualAdjust :: Sequence a -> Sequence a-qualAdjust (Seq _ _ Nothing) = error "no quality data - impossible!"-qualAdjust sq@(Seq l d (Just q)) =  Seq l_trim (B.unfoldr conv avgs) (Just q)-    where avgs = (sliding_avg 1 q, sliding_avg 20 q, d)-          conv (a:as,s:ss,dd) = Just (if a>20 && s>25 then toUpper (B.head dd)+--   Upper case is >limit and sliding avg >limit*1.3+--   Trimming suggestions so that ends with average windows q<limit are clipped+qualAdjust :: Double -> Sequence a -> Sequence a+qualAdjust _ (Seq _ _ Nothing) = error "no quality data - impossible!"+qualAdjust limit sss =  Seq l_trim (B.unfoldr conv avgs) (Just q)+    where sq@(Seq l d (Just q)) = seqmap (\(cv,qv) -> (cv,if cv `elem` "ACGTacgt" then qv else 0)) sss+          avgs = (sliding_avg 1 q, sliding_avg 20 q, d)+          conv (a:as,s:ss,dd) = Just (if a>limit && s>limit*1.3 then toUpper (B.head dd)                                       else if  a<4 || s<7 then 'n'                                            else toLower (B.head dd),(as,ss,B.tail dd))           conv ([],[],dd) | B.null dd = Nothing -- else broken invariant           conv _ = error "internal error in 'qualAdjust/conv'" -          (trim_left,trim_right) = let (_,as,_) = avgs -                                   in (length $ takeWhile (<20) as, fromIntegral (seqlength sq) - length (takeWhile (<20) $ reverse as))-          l_trim = B.concat (l:map B.pack [" QTRIM: ",show trim_left," ",show trim_right])+          (trim_left,trim_right) = let (vs,as,_) = avgs+                                       atrim_l = length $ takeWhile (<limit) as+                                       atrim_r = length (takeWhile (<limit) $ reverse as)+                                   in ( atrim_l + length (takeWhile (<(limit*1.3)) (drop atrim_l vs))+                                      , fromIntegral (seqlength sq) - (atrim_r + length (takeWhile (<limit*1.3) $ drop atrim_r $ reverse vs)))+          l_trim = if trim_left > trim_right then B.concat [l,B.pack " QTRIM: 0 0"]+                   else B.concat (l:map B.pack [" QTRIM: ",show trim_left," ",show trim_right])  myReadPhd :: FilePath -> IO (FilePath,Sequence Nuc) myReadPhd f = unsafeInterleaveIO (do p <- readPhd f@@ -242,11 +389,7 @@ plot term (f,s) = case seqlength s of     0 -> hPutStrLn stderr ("cannot plot empty sequence: "++f) >> return ()     _ -> do     (i,o,e,p) <- runInteractiveCommand "gnuplot -persist"-                hPutStr i header-                hPutStr i $ unlines $ map (show . ord) $ toStr $ seqqual s-                hPutStr i "e\n"-                hPutStr i $ unlines $ map show $ sliding_avg 20 $ seqqual s-                hPutStr i "e\n"+                hPutStr i (mkGnuplot term (f,s))                 hClose i                 x <- waitForProcess p                 hGetContents o >>= hPutStr stderr@@ -254,14 +397,24 @@                 case x of ExitSuccess ->  return ()                           ExitFailure j -> hPutStr stderr (errmsg++show j)                                            >> return ()+    where errmsg = "'gnuplot' failed for "++f++" with exit code "++-- | Build a GNUplot file for graphing quality+mkGnuplot :: (FilePath -> String) -> (FilePath,Sequence Nuc) -> String+mkGnuplot term (f,s) = concat +                [ header+                , unlines $ map (show . ord) $ toStr $ seqqual s+                , "e\n"+                , unlines $ map show $ sliding_avg 20 $ seqqual s+                , "e\n"+                ]     where name = toStr (seqlabel s)           arrow (i,c) = "set arrow "++show i++" from "++show c++",0 to "++show c ++",20 nohead\n"           mtrim =  concatMap arrow $ zip [(1::Int)..] (fst $ trims s)           header = term f                    ++"set xlabel \"position\"\nset ylabel \"quality\"\n"                    ++ mtrim-                   ++ "set title \""++name++"\"\nset yrange [0:100]\nplot \"-\" t \"qual\" with points,\"-\" t \"avg\" with lines lc 0, 20 with lines t \"thresh\"\n"-          errmsg = "'gnuplot' failed for "++f++" with exit code "+                   ++ "set title \""++name++"\"\nset yrange [0:100]\nplot \"-\" t \"qual\" with points,\"-\" t \"avg\" with lines lt 3, 20 with lines lt 2 t \"thresh\"\n"  -- | Look for trimming information.  Phred outputs "TRIM:" in the --   sequence header, followed by trimming information.  Lucy outputs@@ -288,12 +441,13 @@ test_seq :: Sequence Nuc                                      test_seq = Seq (B.pack "foo TRIM: 1 10 QTRIM: 4 15") (B.pack "1234567890abcdefghij") Nothing --- | Calculate a sliding average. Slightly biased for even i. sliding_avg :: Int -> QualData -> [Double]-sliding_avg i q = map avg $ take (fromIntegral $ B.length q) (iblocks i qs)-    where avg xs = fromIntegral (sum xs) / fromIntegral (length xs)-          qs = map ord . toStr $ q---- | select centered words-iblocks :: Int -> [a] -> [[a]]-iblocks i ss = [take n ss | n <- [(i+1) `div` 2..(i-1)]] ++ map (take i) (tails ss)+sliding_avg i q = let k = max 1 (min i (fromIntegral $ B.length q) `div` 2)+                      sum1 = fromIntegral $ sum $ map ord $ take k $ toStr $ q+                  in take (fromIntegral $ B.length q) $ unfoldr getAvg (sum1,k,B.drop (fromIntegral k) q,q)+    where getAvg :: (Int,Int,QualData,QualData) -> Maybe (Double,(Int,Int,QualData,QualData))+          getAvg (!a,!n,!q1,!q2) | B.null q1 && B.null q2 = Nothing+                                 | B.null q1 = let x = a-ord (B.head q2) in Just (fromIntegral a/fromIntegral n,(x, n-1, q1, B.tail q2))+                                 | n<i       = let x = a+ord (B.head q1) in Just (fromIntegral a/fromIntegral n,(x, n+1, B.tail q1, q2))+                                 | otherwise = let x = a+ord (B.head q1)-ord (B.head q2) +                                               in Just (fromIntegral a/fromIntegral n, (x, n, B.tail q1, B.tail q2))
+ src/FakeQual.hs view
@@ -0,0 +1,33 @@+{-| This program reads a Fasta file, and generates a corresponding qual+    file with made-up quality values.  Useful when you want to combine +    sequences with and without quality in the same analysis.+-}++{-# Language DeriveDataTypeable #-}++module Main where++import System.Console.CmdArgs+import qualified Data.ByteString.Lazy as B+import System.IO (stdin,stdout)+import Bio.Sequence++version = "fakequal v0.1.4, copyright 2010 Ketil Malde"++data Opts = O { qval :: Int, input :: [FilePath] } +          deriving (Show,Data,Typeable)++modes = mode $ O { qval = 15 &= text "value to use for quality" & typ "Int" & flag "q"+                 , input = def &= args & typFile }+                 &= prog "fakequal" +                 & text "generate fake quality information for fasta files"++main = do+  cf <- cmdArgs version [modes]+  -- print cf+  ss <- if null (input cf) then hReadFasta stdin+        else concat `fmap` mapM readFasta (input cf)+  hWriteQual stdout (map (addqual $ qval cf) ss)++addqual :: Int -> Sequence t -> Sequence t+addqual q (Seq h s _) = Seq h s (Just $ B.map (const $ fromIntegral q) s)