diff --git a/seqloc-datafiles.cabal b/seqloc-datafiles.cabal --- a/seqloc-datafiles.cabal +++ b/seqloc-datafiles.cabal @@ -1,5 +1,5 @@ Name: seqloc-datafiles -Version: 0.1 +Version: 0.1.1 Cabal-Version: >= 1.4 Synopsis: Read and write BED and GTF format genome annotations Description: Read and write BED and GTF format genome annotations @@ -18,7 +18,7 @@ Library Exposed-modules: Bio.SeqLoc.Bed, Bio.SeqLoc.GTF, Bio.SeqLoc.TranscriptTable Other-modules: Bio.SeqLoc.ZeptoUtils - Build-depends: base >= 4.3 && < 5, bytestring, haskell98, + Build-depends: base >= 4.2 && < 5, bytestring, haskell98, attoparsec >= 0.8.5, hashable, unordered-containers, iteratee >= 0.8.1, seqloc >= 0.0 Hs-Source-Dirs: src @@ -27,11 +27,11 @@ Executable gtf-to-bed Main-is: GtfToBed.hs Other-modules: Bio.SeqLoc.GTF, Bio.SeqLoc.Bed, Bio.SeqLoc.ZeptoUtils - Build-depends: base >= 4.3 && < 5, bytestring, haskell98, + Build-depends: base >= 4.2 && < 5, bytestring, haskell98, attoparsec >= 0.8.5, hashable, unordered-containers, iteratee >= 0.8.1, seqloc >= 0.0, transformers, monads-tf Hs-Source-Dirs: src - Ghc-options: -Wall + Ghc-options: -Wall -rtsopts C-Sources: src/rtsopts.c Executable test-gtf @@ -39,7 +39,7 @@ Buildable: False Main-is: TestGtf.hs Other-modules: Bio.SeqLoc.GTF, Bio.SeqLoc.TranscriptTable, Bio.SeqLoc.ZeptoUtils - Build-depends: base >= 4.3 && < 5, bytestring, haskell98, + Build-depends: base >= 4.2 && < 5, bytestring, haskell98, attoparsec >= 0.8.5, hashable, unordered-containers, iteratee >= 0.8.1, seqloc >= 0.0, QuickCheck, random @@ -52,7 +52,7 @@ Buildable: False Main-is: TestBed.hs Other-modules: Bio.SeqLoc.Bed, Bio.SeqLoc.TranscriptTable, Bio.SeqLoc.ZeptoUtils - Build-depends: base >= 4.3 && < 5, bytestring, haskell98, + Build-depends: base >= 4.2 && < 5, bytestring, haskell98, attoparsec >= 0.8.5, hashable, unordered-containers, iteratee >= 0.8.1, seqloc >= 0.0, QuickCheck, random @@ -64,7 +64,7 @@ Buildable: False Main-is: GtfIntrons.hs Other-modules: Bio.SeqLoc.GTF, Bio.SeqLoc.TranscriptTable, Bio.SeqLoc.ZeptoUtils - Build-depends: base >= 4.3 && < 5, bytestring, haskell98, + Build-depends: base >= 4.2 && < 5, bytestring, haskell98, attoparsec >= 0.8.5, hashable, unordered-containers, iteratee >= 0.8.1, seqloc >= 0.0, QuickCheck, random diff --git a/src/Bio/SeqLoc/GTF.hs b/src/Bio/SeqLoc/GTF.hs --- a/src/Bio/SeqLoc/GTF.hs +++ b/src/Bio/SeqLoc/GTF.hs @@ -3,7 +3,7 @@ {-| Utilities for reading and writing GTF format gene annotations -} module Bio.SeqLoc.GTF - ( readGtfTranscripts + ( readGtfTranscripts, readGtfTranscriptsErr , transcriptToGtf ) where @@ -73,6 +73,12 @@ where gtfTrxsIter = Iter.joinI . IterChar.enumLinesBS . Iter.joinI . gtflineIter $ Iter.foldl' insertGtfLine trxs0 trxs0 = GtfTrxs HM.empty HM.empty HM.empty +readGtfTranscriptsErr :: FilePath -> IO ([Transcript], [String]) +readGtfTranscriptsErr = Iter.fileDriver gtfTrxsIter >=> + return . mkTranscriptsErr + where gtfTrxsIter = Iter.joinI . IterChar.enumLinesBS . Iter.joinI . gtflineIter $ Iter.foldl' insertGtfLine trxs0 + trxs0 = GtfTrxs HM.empty HM.empty HM.empty + mkTranscripts :: GtfTrxs -> Either String [Transcript] mkTranscripts trxs = go [] allTrxs where allTrxs = HM.toList . gtfTogene $ trxs @@ -84,10 +90,23 @@ where exons = fromMaybe [] . HM.lookup trxname . gtfExonLocs $ trxs cdses = fromMaybe [] . HM.lookup trxname . gtfCdsLocs $ trxs +mkTranscriptsErr :: GtfTrxs -> ([Transcript], [String]) +mkTranscriptsErr trxs = go ([], []) allTrxs + where allTrxs = HM.toList . gtfTogene $ trxs + go curr [] = curr + go (currtrx, currerr) (trxAndGene:rest) = case mkOne trxAndGene of + Left err -> let nexterr = err : currerr + in nexterr `seq` go (currtrx, nexterr) rest + Right t -> let next = t : currtrx + in next `seq` go (next, currerr) rest + mkOne (trxname, genename) = mkTranscript trxname exons cdses genename + where exons = fromMaybe [] . HM.lookup trxname . gtfExonLocs $ trxs + cdses = fromMaybe [] . HM.lookup trxname . gtfCdsLocs $ trxs + mkTranscript :: BS.ByteString -> [ContigSeqLoc] -> [ContigSeqLoc] -> BS.ByteString -> Either String Transcript mkTranscript trx exons cdses gene = moderr $ do loc <- exonLoc exons cdsloc <- cdsLoc loc cdses - return $! Transcript (SeqName gene) (SeqName trx) loc cdsloc + return $ Transcript (SeqName gene) (SeqName trx) loc cdsloc where moderr = either (Left . (("Transcript " ++ show trx ++ ": ") ++)) Right exonLoc :: [ContigSeqLoc] -> Either String SpliceSeqLoc diff --git a/src/GtfToBed.hs b/src/GtfToBed.hs --- a/src/GtfToBed.hs +++ b/src/GtfToBed.hs @@ -12,6 +12,10 @@ import Bio.SeqLoc.Bed import Bio.SeqLoc.GTF +import qualified Bio.SeqLoc.Location as Loc +import Bio.SeqLoc.LocRepr +import Bio.SeqLoc.Strand +import Bio.SeqLoc.Transcript main :: IO () main = getArgs >>= handleOpt . getOpt RequireOrder optDescrs @@ -23,29 +27,51 @@ hPutStrLn stderr errs doGtfToBed :: FilePath -> Conf -> IO () -doGtfToBed gtf conf = readGtfTranscripts gtf >>= writeTranscriptsOut +doGtfToBed gtf conf = readGtfTranscriptsErr gtf >>= handleErrors >>= writeTranscriptsOut . map adjustStop where writeTranscriptsOut trxs = withOutHandle conf $ \h -> mapM_ (BS.hPutStrLn h . transcriptToBedStd) trxs + adjustStop = if confStopIncluded conf then removeDoubleStop else id + removeDoubleStop trx = trx { cds = cds trx >>= shortenLoc } + shortenLoc l = Loc.clocOutof (Loc.fromBoundsStrand 0 end Fwd) l + where end = max 0 $ Loc.length l - 4 + handleErrors (trxs, errs) = do maybe putErrs filePutErrs $ confBadTranscripts conf + return trxs + where putErrs = hPutErrs stderr + filePutErrs badtx = withFile badtx WriteMode hPutErrs + hPutErrs h = hPutStr h (unlines errs) withOutHandle :: Conf -> (Handle -> IO a) -> IO a withOutHandle conf m = maybe (m stdout) (\outname -> withFile outname WriteMode m) $ confOutput conf data Conf = Conf { confOutput :: !(Maybe FilePath) + , confBadTranscripts :: !(Maybe FilePath) + , confStopIncluded :: !Bool } deriving (Show) data Arg = ArgOutput { unArgOutput :: !String } + | ArgBadTranscripts { unArgBadTranscripts :: !String } + | ArgStopIncluded deriving (Show, Read, Eq, Ord) argOutput :: Arg -> Maybe String argOutput (ArgOutput del) = Just del argOutput _ = Nothing +argBadTranscripts :: Arg -> Maybe String +argBadTranscripts (ArgBadTranscripts bad) = Just bad +argBadTranscripts _ = Nothing + optDescrs :: [OptDescr Arg] -optDescrs = [ Option ['o'] ["output"] (ReqArg ArgOutput "OUTFILE") "Output filename" +optDescrs = [ Option ['o'] ["output"] (ReqArg ArgOutput "OUTFILE") "Output filename" + , Option ['b'] ["bad-transcripts"] (ReqArg ArgBadTranscripts "BADFILE") "Write bad transcript list to file" + , Option ['s'] ["stop-included"] (NoArg ArgStopIncluded) "GTF file includes stop codon in CDS" ] argsToConf :: [Arg] -> Either String Conf argsToConf = runReaderT conf where conf = Conf <$> - findOutput + findOutput <*> + findBadTranscripts <*> + (ReaderT $ return . elem ArgStopIncluded) findOutput = ReaderT $ return . listToMaybe . mapMaybe argOutput + findBadTranscripts = ReaderT $ return . listToMaybe . mapMaybe argBadTranscripts