diff --git a/dephd.cabal b/dephd.cabal --- a/dephd.cabal +++ b/dephd.cabal @@ -1,5 +1,5 @@ Name: dephd -Version: 0.1.1 +Version: 0.1.3 License: GPL License-File: LICENSE @@ -7,14 +7,22 @@ Maintainer: Ketil Malde Category: Bioinformatics -Synopsis: Analyze 'phred' output (.phd files) +Synopsis: Analyze quality of nucleotide sequences. Description: dephd - A simple tool for base calling and quality appraisal. . Reads files in phd-format (phred output), either specified individually, - or in a directory (use the --dir option to read directories). + or in a directory (use the --input-dirs option to read directories + or --input-list to read from an index file). Can also read FASTA with an + associated quality file. + . + 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. + . The Darcs repository is at . -Build-Depends: base>3, bio >= 0.3.3.4, regex-compat, bytestring, process, directory +HomePage: http://malde.org/~ketil/biohaskell/dephd +Build-Depends: base>=3 && <4, bio >= 0.4, regex-compat, bytestring, process, directory Build-Type: Simple Tested-with: GHC==6.8.3 @@ -22,4 +30,4 @@ Executable: dephd Main-Is: Dephd.hs Hs-Source-Dirs: src -Ghc-Options: -Wall -O2 -funbox-strict-fields +Ghc-Options: -Wall diff --git a/src/Dephd.hs b/src/Dephd.hs --- a/src/Dephd.hs +++ b/src/Dephd.hs @@ -28,15 +28,17 @@ -- Option Handling -- ------------------------------------------------------------ -data MyOpts = O { actions :: [(FilePath,Sequence) -> IO ()] -- ^ Apply to each sequence +data MyOpts = O { actions :: [(FilePath,Sequence Nuc) -> IO ()] -- ^ Apply to each sequence , outputs :: [Handle] -- ^ Output handles that must be closed - , filters :: (FilePath,Sequence) -> (FilePath,Sequence) -- ^ Filter sequences before processing - , inputs :: [String] -> IO [(FilePath,Sequence)] -- ^ Turn args into sequences + , filters :: (FilePath,Sequence Nuc) -> (FilePath,Sequence Nuc) -- ^ Filter sequences before processing + , zerofilter :: (FilePath,Sequence Nuc) -> Bool + , inputs :: [String] -> IO [(FilePath,Sequence Nuc)] -- ^ Turn args into sequences , verbose :: Bool -- ^ Verbose output (progress reporting) and sequence trimming } defaultopts :: MyOpts defaultopts = O { actions = [], outputs = [], filters = id + , zerofilter = const True , inputs = readFiles, verbose = False } getOptions :: IO (MyOpts, [String], [String]) @@ -61,13 +63,14 @@ , Option ['X'] ["output-xplot"] (NoArg (setplot X)) "Display quality plots" -- Filter options - , Option ['t'] ["filter-trim"] (NoArg filterTrim) "Trim output sequences" + , 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" -- 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" , Option ['i'] ["input-fasta-qual"] (NoArg inputFQ) "Read a fasta- and optionally a qual file" + , Option ['z'] ["ignore-empty-seqs"] (NoArg (\opt -> return opt { zerofilter = zeroFilter})) "Eliminate zero-length sequences from output" ] where inputList opt = return $ opt { inputs = \args -> case args of @@ -78,10 +81,11 @@ { inputs = \arg -> do ss <- case arg of [fa,q] -> readFastaQual fa q [fa] -> readFasta fa _ -> error ("Too many files specified for fasta/qual input: "++show arg) - return $ map (\s -> (toStr $ seqlabel s, s)) ss } + 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)) + zeroFilter = (/=0). seqlength . snd addAction act arg opt = let as = actions opt hs = outputs opt condOpenArg fn @@ -108,7 +112,7 @@ 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 :: (FilePath,Sequence Nuc) -> (FilePath,Sequence Nuc) 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) @@ -124,7 +128,7 @@ hasgp :: IO Bool hasgp = return . isJust =<< findExecutable "gnuplot" -readDirs, readFiles :: [FilePath] -> IO [(FilePath,Sequence)] +readDirs, readFiles :: [FilePath] -> IO [(FilePath,Sequence Nuc)] readDirs dirs = mapM' myReadPhd =<< filterM isPhdFile =<< return . concat =<< mapM myGetDirectoryContents =<< filterM doesDirectoryExist dirs readFiles = mapM' myReadPhd @@ -148,7 +152,7 @@ let process xs = if not (verbose opts) then return xs else countIO ("processing "++show (length xs)++" sequences: ") ", done.\n" 100 xs - mapM_ (sequence_ . zipWith ($) (actions opts) . repeat) =<< return . map (filters opts) =<< process =<< (inputs opts) fs + mapM_ (sequence_ . zipWith ($) (actions opts) . repeat) =<< return . filter (zerofilter opts) . map (filters opts) =<< process =<< (inputs opts) fs mapM_ hClose (outputs opts) phd_rx :: Regex @@ -164,9 +168,9 @@ -- | Adjust sequence content according to quality. -- Upper case is >20 and sliding avg >25 -qualAdjust :: Sequence -> Sequence +qualAdjust :: Sequence a -> Sequence a qualAdjust (Seq _ _ Nothing) = error "no quality data - impossible!" -qualAdjust (Seq l d (Just q)) = Seq l (B.unfoldr conv avgs) (Just q) +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) else if a<4 || s<7 then 'n' @@ -174,7 +178,11 @@ conv ([],[],dd) | B.null dd = Nothing -- else broken invariant conv _ = error "internal error in 'qualAdjust/conv'" -myReadPhd :: FilePath -> IO (FilePath,Sequence) + (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]) + +myReadPhd :: FilePath -> IO (FilePath,Sequence Nuc) myReadPhd f = unsafeInterleaveIO (do p <- readPhd f p `seq` return (f,p)) @@ -205,7 +213,7 @@ myhead x = if null x then 0 else head x -- | Report (various?) quality estimates -qualchk :: Sequence -> [String] +qualchk :: Sequence a -> [String] qualchk s = [toStr $ seqlabel s ,printf "%.1f" avgqual] ++ map show [s15, s30, a20] ++[show qtot] @@ -226,11 +234,11 @@ sortOn f = sortBy (\x y -> compare (f x) (f y)) -- | Plot the quality of a sequence in the background -bgplot :: (FilePath -> String) -> (FilePath,Sequence) -> IO ThreadId +bgplot :: (FilePath -> String) -> (FilePath,Sequence Nuc) -> IO ThreadId bgplot x = forkIO . plot x -- | Feed the quality data to gnuplot (check if installed?) -plot :: (FilePath -> String) -> (FilePath,Sequence) -> IO () +plot :: (FilePath -> String) -> (FilePath,Sequence Nuc) -> 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" @@ -257,18 +265,28 @@ -- | 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. +-- just a bunch of numbers, the latter two are trimming information. -- Returns the trim points, and the header with trimming info removed. -trims :: Sequence -> ([Int],SeqData) +trims :: Sequence Nuc -> ([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 + (pt,rest1) = get_trim (B.pack "TRIM:") ws + (dt,rest2) = get_trim (B.pack "QTRIM:") rest1 + mintrims [a,b] [c,d] = [max a c,min b d] + mintrims [] x = x + mintrims x _ = x + in 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) + else (mintrims pt dt, B.unwords rest2) + +get_trim :: B.ByteString -> [B.ByteString] -> ([Int],[B.ByteString]) +get_trim key sh = case dropWhile ((/=) key) sh of + (_:x:y:rest) -> ([read $ B.unpack x, read $ B.unpack y] + ,takeWhile ((/=) key) sh ++ rest) + _ -> ([],sh) + +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]