diff --git a/README b/README --- a/README +++ b/README @@ -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 diff --git a/dephd.cabal b/dephd.cabal --- a/dephd.cabal +++ b/dephd.cabal @@ -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 . + The Darcs repository is at . -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 diff --git a/src/Dephd.hs b/src/Dephd.hs --- a/src/Dephd.hs +++ b/src/Dephd.hs @@ -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)