dephd 0.1.1 → 0.1.3
raw patch · 2 files changed
+53/−27 lines, 2 filesdep ~basedep ~bio
Dependency ranges changed: base, bio
Files
- dephd.cabal +13/−5
- src/Dephd.hs +40/−22
dephd.cabal view
@@ -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 <ketil@malde.org> 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 <http://malde.org/~ketil/biohaskell/dephd>. -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
src/Dephd.hs view
@@ -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]