dephd 0.0 → 0.1
raw patch · 3 files changed
+84/−69 lines, 3 filesdep ~bytestringnew-uploader
Dependency ranges changed: bytestring
Files
- README +25/−26
- dephd.cabal +3/−3
- src/Dephd.hs +56/−40
README view
@@ -26,39 +26,38 @@ Usage ----- - dephd --rank files..- dephd --rank --dir dirs..--Outputs (to standard output) a summary of all phd files, including-sequence name, average quality, length of longest contiguous run with-qualities >= 15, 20 and 30, and longest run with sliding average-quality 20 or better.-- dephd --call files..- dephd --call --dir dirs..+A brief usage report is printed if you run 'dephd -h'. Somewhat more detailed: -Produce files 'dephd.fasta' and 'dephd.qual' in the current directory,-containing sequence and quality data, respectively. Bases estimated-(currently very conservatively) to be of good quality are in upper-case, very low quality is output as lower case 'n's.+Input is specified either as a list of phd-files (typcially generated+by Phred), a list of directories containing phd-files (using the+--input-dirs) option, a file containing a list of names of phd-files+(--input-list), or a Fasta and associated quality file (-i foo.fasta+foo.qual). - dephd --plot files..- dephd --plot --dir dirs..- dephd --plot -X files..- dephd --plot -X --dir dirs..+Output is specified by -J, -X, -P, -R foo.ranks, -F foo.fasta, and/or+-Q foo.qual. The first three generate a plot of sequence quality in+JPEG files, an X window, or Postscript files, respectively. If you+use -X on multiple files, hit q to terminate one window and go to the+next. -Produce a plot of sequence quality, along with a sliding average.-With -X display it directory, without -X, produce a jpg file with the-plot. A similar option --plotall generates and displays all plots-directly, instead of one at a time (only useful with -X as well).+The remaining three 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. +Filtering can be specified with the -t option, which interprets+trimming information from Phred or Lucy, and chops off the offending+parts, or with the -q options, which masks poor quality parts of+sequences to lower case, and really poor quality parts to 'n'+characters. Bugs ---- -Not many, I hope. Specifying more than one action at a time will pull-all sequences into memory, but a single action should stream okay.-Approx 15K phd-files can be --call'ed OR --plot'ed OR --rank'ed with-less than 100Mb of RAM.+Not many, I hope. The program should work in (approximately) constant+space, and be able to deal with large amounts of sequences. For further questions, email me at <ketil@malde.org>
dephd.cabal view
@@ -1,5 +1,5 @@ Name: dephd-Version: 0.0+Version: 0.1 License: GPL License-File: LICENSE @@ -12,9 +12,9 @@ . Reads files in phd-format (phred output), either specified individually, or in a directory (use the --dir option to read directories).- The Darcs repository is at <http://malde.org/~ketil/dephd>.+ The Darcs repository is at <http://malde.org/~ketil/biohaskell/dephd>. -Build-Depends: base>3, bio >= 0.3, regex-compat, bytestring==0.9.0.1, process, directory+Build-Depends: base>3, bio >= 0.3, regex-compat, bytestring, process, directory Build-Type: Simple Tested-with: GHC==6.8.3
src/Dephd.hs view
@@ -8,7 +8,7 @@ import Control.Concurrent import Control.Monad import Data.Char-import Data.List (tails,groupBy,sortBy)+import Data.List (tails,groupBy,sortBy,isPrefixOf) import Data.Maybe import System.Console.GetOpt import System.Directory@@ -22,6 +22,7 @@ import qualified Data.ByteString.Lazy.Char8 as B import Bio.Sequence+import Bio.Util (countIO) -- ------------------------------------------------------------ -- Option Handling@@ -44,17 +45,20 @@ os' <- foldl (>>=) (return defaultopts) os return (os',ns,es) +data PlotType = J | P | X+ options :: [OptDescr (MyOpts -> IO MyOpts)] options = [ Option ['v'] [] (NoArg (\opt -> return opt {verbose = True})) "Verbose output" , Option ['h'] ["help"] (NoArg (\_ -> do {putStrLn (usage []); exitWith ExitSuccess})) "Display usage information" -- Output options- , Option ['R'] ["output-ranks"] (ReqArg setrank "file") "Set ranked output file"+ , 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 ['P'] ["output-plot"] (NoArg (setplot False)) "Generate quality plots"- , Option ['X'] ["output-xplot"] (NoArg (setplot True)) "Display quality plots"+ , Option ['Q'] ["output-qual"] (ReqArg setqual "file") "Set quality output file"+ , 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" -- Filter options , Option ['t'] ["filter-trim"] (NoArg filterTrim) "Trim output sequences"@@ -80,26 +84,42 @@ filterQual = addFilter (\(p,s) -> (p,qualAdjust s)) addAction act arg opt = let as = actions opt hs = outputs opt- in do h <- openFile arg WriteMode+ condOpenArg fn + | fn == "-" = return stdout + | "-" `isPrefixOf` fn = error ("Refusing output file name ('"++fn+ ++ "') starting with '-'.\n"+ ++ "Use ' ./"++fn++"' if this is what you want.")+ | otherwise = openFile fn WriteMode+ in do h <- condOpenArg arg return opt { actions = act h:as , outputs = h:hs } setrank = addAction (\h -> hPutStrLn h . unwords {- columns -} . qualchk . snd) setfasta = addAction (\h -> hWriteFasta h . return . snd) setqual = addAction (\h -> hWriteQual h . return . snd)- setplot b opt = do++ setplot c opt = do g <- hasgp- if g then let as = actions opt in return $ opt { actions = plot b : as }+ let how = case c of J -> \f -> "set terminal jpeg\nset output \""+ ++subRegex phd_rx f "" ++ ".jpg\"\n"+ P -> \f -> "set terminal postscript color eps\nset output\""+ ++subRegex phd_rx f "" ++ ".eps\"\n"+ X -> const ""+ 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") trimSeq :: (FilePath,Sequence) -> (FilePath,Sequence)-trimSeq (i,s@(Seq h d mq)) = case trims s of- [t1,t2] -> let clip = B.take (fromIntegral t2-fromIntegral t1) . B.drop (fromIntegral t1)- in (i,Seq h (clip d) (case mq of Nothing -> Nothing- Just q -> Just (clip q)))- _ -> (i,s)+trimSeq (i,s@(Seq _ d mq)) = + case trims s of + ([t1,t2],h') -> let clip = B.take (fromIntegral t2-fromIntegral t1) . B.drop (fromIntegral t1)+ s' = Seq h' (clip d) (case mq of Nothing -> Nothing+ Just q -> Just (clip q))+ in (i,appendHeader s' $ unwords ["clipped:",show t1,show t2])+ _ -> (i,s) -- todo: clip only 'n's?+-- todo: add trimming info in header? +-- (Currently, nothing protects agains re-trimming with the same parameters...) hasgp :: IO Bool hasgp = return . isJust =<< findExecutable "gnuplot"@@ -110,7 +130,7 @@ readFiles = mapM' myReadPhd mapM' :: (a -> IO b) -> [a] -> IO [b]-mapM' _ [] = hPutStrLn stderr "Warning: no files found" >> return []+mapM' _ [] = hPutStrLn stderr "Warning: nothing to do!\n(Use '-h' for help)" >> return [] mapM' action xs = mapM action xs -- ------------------------------------------------------------@@ -143,13 +163,13 @@ myGetDirectoryContents d = return . map ((d++"/")++) =<< getDirectoryContents d -- | Adjust sequence content according to quality.--- Upper case is >20 and sliding avg >30+-- Upper case is >20 and sliding avg >25 qualAdjust :: Sequence -> Sequence qualAdjust (Seq _ _ Nothing) = error "no quality data - impossible!" qualAdjust (Seq l d (Just q)) = Seq l (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)- else if a<4 || s<9 then 'n'+ 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'"@@ -206,12 +226,12 @@ sortOn f = sortBy (\x y -> compare (f x) (f y)) -- | Plot the quality of a sequence in the background-bgplot :: Bool -> (FilePath,Sequence) -> IO ThreadId+bgplot :: (FilePath -> String) -> (FilePath,Sequence) -> IO ThreadId bgplot x = forkIO . plot x -- | Feed the quality data to gnuplot (check if installed?)-plot :: Bool -> (FilePath,Sequence) -> IO ()-plot z (f,s) = case seqlength s of+plot :: (FilePath -> String) -> (FilePath,Sequence) -> IO ()+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@@ -227,19 +247,28 @@ ExitFailure j -> hPutStr stderr (errmsg++show j) >> return () where name = toStr (seqlabel s)- f' = subRegex phd_rx f ".jpg" arrow (i,c) = "set arrow "++show i++" from "++show c++",0 to "++show c ++",20 nohead\n"- mtrim = concatMap arrow $ zip [(1::Int)..] (trims s)- header = (if z then "" else "set term jpeg\nset out \""++f'++"\"\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:75]\nplot \"-\" t \"qual\" with points,\"-\" t \"avg\" with lines, 20 with lines t \"thresh\"\n"+ ++ "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 " -trims :: Sequence -> [Int]-trims s = case dropWhile ((/=) (B.pack "TRIM:")) (B.words $ seqheader s) of- (_:x:y:_) -> [read $ B.unpack x, read $ B.unpack y]- _ -> []+-- | Look for trimming information. Phred outputs "TRIM:" in the+-- sequence header, followed by trimming information. Lucy outputs+-- just a bunch of number, the latter two are trimming information.+-- Returns the trim points, and the header with trimming info removed.+trims :: Sequence -> ([Int],SeqData)+trims s = let ws = B.words $ seqheader s+ phred_trim = B.pack "TRIM:"+ in case dropWhile ((/=) phred_trim) ws of+ (_:x:y:rest) -> ([read $ B.unpack x, read $ B.unpack y]+ ,B.unwords (takeWhile ((/=) phred_trim) ws ++ rest))+ _ -> if all (all isDigit . B.unpack) (tail ws) && length ws > 2+ then (map (read . B.unpack) $ reverse $ take 2 $ reverse ws+ ,B.unwords $ reverse $ drop 2 $ reverse ws)+ else ([],seqheader s) -- | Calculate a sliding average. Slightly biased for even i. sliding_avg :: Int -> QualData -> [Double]@@ -250,16 +279,3 @@ -- | 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)---- copied from RBR-countIO :: String -> String -> Int -> [a] -> IO [a]-countIO msg post step xs = sequence $ map unsafeInterleaveIO ((blank >> outmsg (0::Int) >> c):cs)- where (c:cs) = ct 0 xs- output = hPutStr stderr- blank = output ('\r':take 70 (repeat ' '))- outmsg x = output ('\r':msg++show x) >> hFlush stderr- ct s ys = let (a,b) = splitAt (step-1) ys- next = s+step- in case b of [b1] -> map return a ++ [outmsg (s+step) >> hPutStr stderr post >> return b1]- [] -> map return (init a) ++ [outmsg (s+length a) >> hPutStr stderr post >> return (last a)]- _ -> map return a ++ [outmsg s >> return (head b)] ++ ct next (tail b)