diff --git a/README b/README --- a/README +++ b/README @@ -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 diff --git a/dephd.cabal b/dephd.cabal --- a/dephd.cabal +++ b/dephd.cabal @@ -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 . . + Also provides 'fakequal', a utility to generate bogus quality values, + which are sometimes needed by less flexible tools. + . The Darcs repository is at . 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 diff --git a/src/Dephd.hs b/src/Dephd.hs --- a/src/Dephd.hs +++ b/src/Dephd.hs @@ -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 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 ( 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 Sequence t -> Sequence t +addqual q (Seq h s _) = Seq h s (Just $ B.map (const $ fromIntegral q) s)