diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License + +Copyright (c) 2008-2011 Nicholas Ingolia + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/seqloc-datafiles.cabal b/seqloc-datafiles.cabal new file mode 100644 --- /dev/null +++ b/seqloc-datafiles.cabal @@ -0,0 +1,73 @@ +Name: seqloc-datafiles +Version: 0.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 +License: MIT +License-File: LICENSE +Author: Nick Ingolia +Maintainer: nick@ingolia.org +Build-Type: Simple + +Category: Bioinformatics + +Flag Tests + Description: Build test program + Default: False + +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, + attoparsec >= 0.8.5, hashable, unordered-containers, + iteratee >= 0.8.1, seqloc >= 0.0 + Hs-Source-Dirs: src + Ghc-options: -Wall + +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, + attoparsec >= 0.8.5, hashable, unordered-containers, + iteratee >= 0.8.1, seqloc >= 0.0, transformers, monads-tf + Hs-Source-Dirs: src + Ghc-options: -Wall + C-Sources: src/rtsopts.c + +Executable test-gtf + if !flag(Tests) + 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, + attoparsec >= 0.8.5, hashable, unordered-containers, + iteratee >= 0.8.1, seqloc >= 0.0, + QuickCheck, random + Hs-Source-Dirs: src, test + Ghc-options: -Wall -rtsopts + C-Sources: src/rtsopts.c + +Executable test-bed + if !flag(Tests) + 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, + attoparsec >= 0.8.5, hashable, unordered-containers, + iteratee >= 0.8.1, seqloc >= 0.0, + QuickCheck, random + Hs-Source-Dirs: src, test + Ghc-options: -Wall + +Executable gtf-introns + if !flag(Tests) + 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, + attoparsec >= 0.8.5, hashable, unordered-containers, + iteratee >= 0.8.1, seqloc >= 0.0, + QuickCheck, random + Hs-Source-Dirs: src, test + Ghc-options: -Wall -rtsopts + C-Sources: src/rtsopts.c diff --git a/src/Bio/SeqLoc/Bed.hs b/src/Bio/SeqLoc/Bed.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/Bed.hs @@ -0,0 +1,144 @@ +{-# LANGUAGE OverloadedStrings #-} + +{-| Utilities for reading and writing BED format gene annotations -} +module Bio.SeqLoc.Bed + ( readBedTranscripts + , bedZP, bedTranscriptEnum + , transcriptToBed, transcriptToBedStd + ) + where + +import Control.Applicative +import Control.Monad +import qualified Data.ByteString.Char8 as BS +import Data.List +import Data.Ord + +import qualified Data.Attoparsec.Zepto as ZP +import qualified Data.Iteratee as Iter +import qualified Data.Iteratee.Char as IterChar + +import Bio.SeqLoc.LocRepr +import qualified Bio.SeqLoc.Location as Loc +import Bio.SeqLoc.OnSeq +import qualified Bio.SeqLoc.Position as Pos +import qualified Bio.SeqLoc.SpliceLocation as SpLoc +import Bio.SeqLoc.Strand +import Bio.SeqLoc.Transcript + +import Bio.SeqLoc.ZeptoUtils + +-- | Convert a 'Transcript' to a BED annotation line. +transcriptToBedStd :: Transcript -> BS.ByteString +transcriptToBedStd = transcriptToBed "0" "0" + +-- | Convert a 'Transcript' to a BED annotation line, specifying the +-- /score/ and /itemRGB/ fields. +transcriptToBed :: BS.ByteString -- ^ score + -> BS.ByteString -- ^ itemRGB + -> Transcript -- ^ transcript + -> BS.ByteString +transcriptToBed score rgb trx = unfields fields + where unfields = BS.intercalate (BS.singleton '\t') + fields = [ unSeqName chrom + , repr $ chromStart + , repr $ chromEnd + 1 + , unSeqName . trxId $ trx + , score + , strandchr + , repr $ thickStart + , repr $ thickEnd + 1 + , rgb + , BS.pack . show . length $ blockSizes + , unCommaList blockSizes + , unCommaList blockStarts + ] + (OnSeq chrom loc) = location trx + (chromStart, chromEnd) = Loc.bounds loc + strandchr = case Loc.strand loc of Fwd -> "+"; RevCompl -> "-" + (thickStart, thickEnd) = maybe noCds (Loc.bounds . unOnSeq) . cdsLocation $ trx + noCds = (chromStart, chromStart - 1) + contigs = sortBy (comparing Loc.offset5) . Loc.toContigs $ loc + blockSizes = map Loc.length contigs + blockStarts = map (subtract chromStart . Loc.offset5) contigs + unCommaList = BS.concat . map (flip BS.append (BS.singleton ',') . repr) + +-- | Read all BED format annotations in a BED file +readBedTranscripts :: FilePath -> IO [Transcript] +readBedTranscripts = Iter.fileDriver (bedTranscriptEnum Iter.stream2list) + +-- | Iteratee to convert an 'Iter.Iteratee' over a 'BS.ByteString', +-- such as the standard 'Iter.fileDriver', into an iteratee over a +-- list of 'Transcript' annotations from the file. +bedTranscriptEnum :: (Monad m) => Iter.Iteratee [Transcript] m a -> Iter.Iteratee BS.ByteString m a +bedTranscriptEnum = Iter.joinI . IterChar.enumLinesBS . Iter.joinI . bedLineEnum + +bedLineEnum :: (Monad m) => Iter.Enumeratee [BS.ByteString] [Transcript] m a +bedLineEnum = Iter.convStream $ Iter.head >>= liftM (: []) . handleErr . ZP.parse bedZP + where handleErr = either (Iter.throwErr . Iter.iterStrExc) return + +-- | Minimalistic 'ZP.Parser'-style parser for a BED format line, not +-- including the trailing newline. +bedZP :: ZP.Parser Transcript +bedZP = do chrom <- field -- The name of the chromosome + chromStart <- decfield -- The starting position of the + -- feature in the chromosome or + -- scaffold. The first base in a + -- chromosome is numbered 0 + chromEnd <- decfield -- The ending position of the feature + -- in the chromosome or scaffold. The + -- chromEnd base is not included in + -- the display of the feature. For + -- example, the first 100 bases of a + -- chromosome are defined as + -- chromStart=0, chromEnd=100, and + -- span the bases numbered 0-99. + name <- field -- Defines the name of the BED line. + _score <- dropField -- A score between 0 and 1000. + str <- strand -- Defines the strand + thickStart <- decfield -- The starting position at which + -- the feature is drawn thickly (for + -- example, the start codon in gene + -- displays). + thickEnd <- decfield -- The ending position at which the + -- feature is drawn thickly (for + -- example, the stop codon in gene + -- displays). + _itemRGB <- dropField -- An RGB value of the form R,G,B + -- (e.g. 255,0,0). + blockCount <- decfield -- The number of blocks (exons) in + -- the BED line. + blockSizes <- commas blockCount decimal <* dropField + -- A comma-separated list of the block sizes. + blockStarts <- commas blockCount decimal + -- A comma-separated list of block starts. + loc <- bedTrxLoc chromStart chromEnd str $ zip blockSizes blockStarts + unless (Loc.bounds loc == (chromStart, chromEnd - 1)) $ + fail $ "Bio.SeqLoc.Bed: bad sploc:" ++ + (BS.unpack . BS.unwords $ [ repr loc, repr chromStart, repr chromEnd ]) + cdsloc <- if thickStart >= thickEnd + then return Nothing + else liftM Just $! bedCdsLoc loc thickStart thickEnd + let n = SeqName $ BS.copy name + c = SeqName $ BS.copy chrom + return $! Transcript n n (OnSeq c loc) cdsloc + +bedTrxLoc :: (Monad m) => Pos.Offset -> Pos.Offset -> Strand -> [(Pos.Offset, Pos.Offset)] -> m SpLoc.SpliceLoc +bedTrxLoc chromStart chromEnd str = maybe badContigs (return . stranded str) . + SpLoc.fromContigs . map blockContig + where blockContig (bsize, bstart) = Loc.fromPosLen (Pos.Pos (chromStart + bstart) Fwd) bsize + badContigs = fail $ "Bio.SeqLoc.Bed: bad blocks in " ++ show (chromStart, chromEnd) + +bedCdsLoc :: (Monad m) => SpLoc.SpliceLoc -> Pos.Offset -> Pos.Offset -> m Loc.ContigLoc +bedCdsLoc loc thickStart thickEnd + = maybe badCdsLoc return $ do + relstart <- Loc.posInto (Pos.Pos thickStart Fwd) loc + relend <- Loc.posInto (Pos.Pos (thickEnd - 1) Fwd) loc + return $! stranded (Loc.strand loc) $ Loc.fromStartEnd (Pos.offset relstart) (Pos.offset relend) + where badCdsLoc = fail $ "Bio.SeqLoc.Bed: bad cds in " ++ + (BS.unpack . BS.unwords $ [ repr loc, repr thickStart, repr thickEnd ]) + +commas :: Int -> ZP.Parser a -> ZP.Parser [a] +commas n p | n < 1 = return [] + | otherwise = (:) <$> p <*> + replicateM (n - 1) (ZP.string "," *> p) diff --git a/src/Bio/SeqLoc/GTF.hs b/src/Bio/SeqLoc/GTF.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/GTF.hs @@ -0,0 +1,187 @@ +{-# LANGUAGE BangPatterns, MagicHash, OverloadedStrings, UnboxedTuples #-} + +{-| Utilities for reading and writing GTF format gene annotations -} + +module Bio.SeqLoc.GTF + ( readGtfTranscripts + , transcriptToGtf + ) + where + +import Control.Applicative +import Control.Monad +import qualified Data.ByteString.Char8 as BS +import Data.ByteString.Internal (c2w) +import Data.List +import Data.Maybe + +import qualified Data.Attoparsec.Char8 as AP (isSpace_w8) +import qualified Data.Attoparsec.Zepto as ZP +import qualified Data.HashMap.Strict as HM +import qualified Data.Iteratee as Iter +import qualified Data.Iteratee.Char as IterChar + +import Bio.SeqLoc.LocRepr +import qualified Bio.SeqLoc.Location as Loc +import Bio.SeqLoc.OnSeq +import qualified Bio.SeqLoc.Position as Pos +import qualified Bio.SeqLoc.SpliceLocation as SpLoc +import Bio.SeqLoc.Strand +import Bio.SeqLoc.Transcript + +import Bio.SeqLoc.ZeptoUtils + +-- | Convert a 'Transcript' to a string consisting of GTF lines. These +-- lines will contain @exon@ lines for the transcript, as well as +-- @CDS@ lines if the 'Transcript' has a 'cds'. +transcriptToGtf :: BS.ByteString -> Transcript -> BS.ByteString +transcriptToGtf src trx = BS.unlines $ exonLines ++ cdsLines + where exonLines = splocLines src trx "exon" (location trx) + cdsLines = maybe [] (splocLines src trx "CDS") $! cdsLocation trx + +splocLines :: BS.ByteString -> Transcript -> BS.ByteString -> SpliceSeqLoc -> [BS.ByteString] +splocLines src trx ftype (OnSeq refname sploc) = map contigLines . Loc.toContigs $ sploc + where contigLines contig = let (start0, end0) = Loc.bounds contig + strchr = case Loc.strand contig of + Fwd -> "+" + RevCompl -> "-" + in unfields [ unSeqName refname + , src + , ftype + , BS.pack . show . Pos.unOffset . (+ 1) $ start0 + , BS.pack . show . Pos.unOffset . (+ 1) $ end0 + , "0.0" + , strchr + , "." + , attrs + ] + attrs = BS.concat [ "gene_id \"", unSeqName . geneId $ trx + , "\"; transcript_id \"", unSeqName . trxId $ trx + , "\"; " ] + unfields = BS.intercalate (BS.singleton '\t') + +-- | Read a GTF annotation file. The entire file is read at once, +-- because a single annotated transcript can span many lines in a GTF +-- file that are not required to occur in any specific order. The +-- transcript 'SpliceSeqLoc' transcript location is assembled from +-- @exon@ annotations and any CDS location is then produced from @CDS@ +-- annotations, with an error occurring if the CDS is not a single +-- contiguous location within the transcript. +readGtfTranscripts :: FilePath -> IO [Transcript] +readGtfTranscripts = Iter.fileDriver gtfTrxsIter >=> + either (ioError . userError) return . mkTranscripts + 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 + go curr [] = Right curr + go curr (trxAndGene:rest) = case mkOne trxAndGene of + Left err -> Left err + Right t -> let next = t : curr in next `seq` go next 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 + where moderr = either (Left . (("Transcript " ++ show trx ++ ": ") ++)) Right + +exonLoc :: [ContigSeqLoc] -> Either String SpliceSeqLoc +exonLoc exons = do + (seqname, rawcontigs) <- allSameName exons + contigs <- sortclocs rawcontigs + sploc <- maybe (badClocs contigs) Right $! SpLoc.fromContigs contigs + return $! OnSeq seqname sploc + where sortclocs clocs = maybe (badClocs clocs) Right . sortContigs $ clocs + badClocs clocs = Left . unwords $ [ "bad transcript exons" ] ++ map (BS.unpack . repr) clocs + +cdsLoc :: SpliceSeqLoc -> [ContigSeqLoc] -> Either String (Maybe Loc.ContigLoc) +cdsLoc _ [] = return Nothing +cdsLoc (OnSeq trxname trxloc) cdses@(_:_) = do + (seqname, rawcontigs) <- allSameName cdses + when (trxname /= seqname) $ Left . unwords $ [ "CDS sequence name mismatch", show trxname, show seqname ] + contigs <- sortclocs rawcontigs + (contigIn0:contigInRest) <- mapM (cdsIntoTranscript trxloc) contigs + cloc <- foldM mergeCLocs contigIn0 contigInRest + return $! Just $! Loc.extend (0, 3) cloc -- Include the stop codon + where sortclocs clocs = maybe badClocs Right . sortContigs $ clocs + where badClocs = Left . unwords $ [ "bad transcript CDSes" ] ++ map (BS.unpack . repr) clocs + +cdsIntoTranscript :: SpLoc.SpliceLoc -> Loc.ContigLoc -> Either String Loc.ContigLoc +cdsIntoTranscript trxloc cdscontig = maybe badIn Right . Loc.clocInto cdscontig $ trxloc + where badIn = Left . unwords $ [ "Mapping CDS contig into transcript: " + , BS.unpack . repr $ cdscontig + , BS.unpack . repr $ trxloc + ] + +mergeCLocs :: Loc.ContigLoc -> Loc.ContigLoc -> Either String Loc.ContigLoc +mergeCLocs cloc0 clocnext + | (Loc.strand cloc0 == Fwd) && (Loc.startPos clocnext == Loc.endPos cloc0 `Pos.slide` 1) + = return $! Loc.extend (0, Loc.length clocnext) cloc0 + | otherwise = Left . unwords $ [ "Merging non-adjacent contigs: " + , BS.unpack . repr $ cloc0 + , BS.unpack . repr $ clocnext + ] + + +allSameName :: [OnSeq a] -> Either String (SeqName, [a]) +allSameName s = case group . map onSeqName $ s of + [(name0:_)] -> return (name0, map unOnSeq s) + names -> Left $ "allSameName: names " ++ show (map (unSeqName . head) names) + +data GtfLine = GtfLine { gtfGene, gtfTrx, gtfFtype :: !BS.ByteString, gtfLoc :: !ContigSeqLoc } deriving (Show) + +data GtfTrxs = GtfTrxs { gtfExonLocs, gtfCdsLocs :: !(HM.HashMap BS.ByteString [ContigSeqLoc]) + , gtfTogene :: !(HM.HashMap BS.ByteString BS.ByteString) + } deriving (Show) + +ftypeCds :: BS.ByteString +ftypeCds = BS.pack "CDS" + +ftypeExon :: BS.ByteString +ftypeExon = BS.pack "exon" + +insertGtfLine :: GtfTrxs -> GtfLine -> GtfTrxs +insertGtfLine trxsin l + | gtfFtype l == ftypeCds = insertCds . insertGene $ trxsin + | gtfFtype l == ftypeExon = insertExon . insertGene $ trxsin + | otherwise = trxsin + where trx = BS.copy $ gtfTrx l + insertGene t0 = {-# SCC "insertGene" #-} t0 { gtfTogene = HM.insert trx (BS.copy . gtfGene $ l) (gtfTogene t0) } + insertCds t0 = {-# SCC "insertCds" #-} t0 { gtfCdsLocs = HM.insertWith (++) trx [gtfLoc l] (gtfCdsLocs t0) } + insertExon t0 = {-# SCC "insertExon" #-} t0 { gtfExonLocs = HM.insertWith (++) trx [gtfLoc l] (gtfExonLocs t0) } + +gtflineIter :: (Monad m) => Iter.Enumeratee [BS.ByteString] [GtfLine] m a +gtflineIter = Iter.convStream $ Iter.head >>= liftM (: []) . handleErr . ZP.parse gtfline + where handleErr = either (Iter.throwErr . Iter.iterStrExc) return + +-- Does NOT consume the remainder of the line +gtfline :: ZP.Parser GtfLine +gtfline = do seqname <- field + _source <- dropField + ftype <- field + start <- decfield + end <- decfield + _score <- dropField + str <- strand + _frame <- dropField + gene <- attr "gene_id" + trx <- attr "transcript_id" + let name = SeqName . BS.copy $ seqname + loc = Loc.fromBoundsStrand (start - 1) (end - 1) str + return $! GtfLine gene trx ftype (OnSeq name loc) + +attr :: String -> ZP.Parser BS.ByteString +attr name = ZP.takeWhile AP.isSpace_w8 *> + ZP.string (BS.pack name) *> + ZP.takeWhile AP.isSpace_w8 *> + ZP.string "\"" *> + ZP.takeWhile (/= c2w '\"') <* + ZP.string "\";" + + + diff --git a/src/Bio/SeqLoc/TranscriptTable.hs b/src/Bio/SeqLoc/TranscriptTable.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/TranscriptTable.hs @@ -0,0 +1,48 @@ +{-| Minimal tab-delimited annotation of 'Transcript' locations. Each +'Transcript' has one line, containing the 'geneId' and 'trxId' fields, +followed by the 'LocRepr' representation of the 'SpliceSeqLoc' +location of the transcript, and then the location of the CDS within +the transcript or "n/a" if there is no CDS. -} + +module Bio.SeqLoc.TranscriptTable + ( readTable, parseLine, writeTable, unparseLine ) + where + +import Control.Monad +import qualified Data.ByteString.Char8 as BS + +import Bio.SeqLoc.LocRepr +import Bio.SeqLoc.OnSeq +import Bio.SeqLoc.Transcript + +-- | Read a transcript table file as a list of annotated transcripts +readTable :: FilePath -> IO [Transcript] +readTable = BS.readFile >=> mapM parseLineErr . BS.lines + where parseLineErr l = maybe err return $! parseLine l + where err = ioError . userError $ "Bad transcript table line " ++ show l + +-- | Parse a transcript table line, not including the newline +parseLine :: BS.ByteString -> Maybe Transcript +parseLine l = case BS.split '\t' l of + [ geneid, trxid, trxlocstr, cdslocstr ] -> do + trxloc <- unreprMaybe trxlocstr + cdsloc <- if cdslocstr == na + then return Nothing + else liftM Just $! unreprMaybe cdslocstr + return $! Transcript (SeqName . BS.copy $ geneid) (SeqName . BS.copy $ trxid) trxloc cdsloc + _ -> Nothing + +-- | Write a transcript table file +writeTable :: FilePath -> [Transcript] -> IO () +writeTable outname = BS.writeFile outname . BS.unlines . map unparseLine + +-- | Produce a single transcript table line for a 'Transcript', newline not included. +unparseLine :: Transcript -> BS.ByteString +unparseLine trx = BS.intercalate (BS.singleton '\t') $ [ unSeqName . geneId $ trx + , unSeqName . trxId $ trx + , repr . location $ trx + , maybe na repr . cds $ trx + ] + +na :: BS.ByteString +na = BS.pack "n/a" diff --git a/src/Bio/SeqLoc/ZeptoUtils.hs b/src/Bio/SeqLoc/ZeptoUtils.hs new file mode 100644 --- /dev/null +++ b/src/Bio/SeqLoc/ZeptoUtils.hs @@ -0,0 +1,32 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Bio.SeqLoc.ZeptoUtils + where + +import Control.Applicative +import qualified Data.ByteString as BSW +import qualified Data.ByteString.Char8 as BS +import Data.ByteString.Internal (c2w) + +import qualified Data.Attoparsec.Char8 as AP (isDigit_w8) +import qualified Data.Attoparsec.Zepto as ZP + +import Bio.SeqLoc.Strand + +strand :: ZP.Parser Strand +strand = ((ZP.string "+\t" *> return Fwd) <|> + (ZP.string "-\t" *> return RevCompl)) + +decfield :: (Integral a) => ZP.Parser a +decfield = decimal <* ZP.string "\t" + +field :: ZP.Parser BS.ByteString +field = ZP.takeWhile (/= c2w '\t') <* ZP.string "\t" + +dropField :: ZP.Parser () +dropField = ZP.takeWhile (/= c2w '\t') *> ZP.string "\t" *> return () + +decimal :: (Integral a) => ZP.Parser a +decimal = decode <$> ZP.takeWhile (AP.isDigit_w8) + where decode = fromIntegral . BSW.foldl' step (0 :: Int) + step a w = a * 10 + fromIntegral (w - 48) diff --git a/src/GtfToBed.hs b/src/GtfToBed.hs new file mode 100644 --- /dev/null +++ b/src/GtfToBed.hs @@ -0,0 +1,51 @@ +module Main + where + +import Control.Applicative +import Control.Monad.Reader +import qualified Data.ByteString.Char8 as BS +import Data.Maybe + +import System.Console.GetOpt +import System.Environment +import System.IO + +import Bio.SeqLoc.Bed +import Bio.SeqLoc.GTF + +main :: IO () +main = getArgs >>= handleOpt . getOpt RequireOrder optDescrs + where handleOpt (_, _, errs@(_:_)) = usage (unlines errs) + handleOpt (args, [gtf], []) = either usage (doGtfToBed gtf) $ argsToConf args + handleOpt (_, _, []) = usage "Specify exactly one GTF file" + usage errs = do prog <- getProgName + hPutStr stderr $ usageInfo prog optDescrs + hPutStrLn stderr errs + +doGtfToBed :: FilePath -> Conf -> IO () +doGtfToBed gtf conf = readGtfTranscripts gtf >>= writeTranscriptsOut + where writeTranscriptsOut trxs = withOutHandle conf $ \h -> + mapM_ (BS.hPutStrLn h . transcriptToBedStd) trxs + +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) + } deriving (Show) + +data Arg = ArgOutput { unArgOutput :: !String } + deriving (Show, Read, Eq, Ord) + +argOutput :: Arg -> Maybe String +argOutput (ArgOutput del) = Just del +argOutput _ = Nothing + +optDescrs :: [OptDescr Arg] +optDescrs = [ Option ['o'] ["output"] (ReqArg ArgOutput "OUTFILE") "Output filename" + ] + +argsToConf :: [Arg] -> Either String Conf +argsToConf = runReaderT conf + where conf = Conf <$> + findOutput + findOutput = ReaderT $ return . listToMaybe . mapMaybe argOutput diff --git a/src/rtsopts.c b/src/rtsopts.c new file mode 100644 --- /dev/null +++ b/src/rtsopts.c @@ -0,0 +1,1 @@ +char *ghc_rts_opts = "-A128m"; diff --git a/test/GtfIntrons.hs b/test/GtfIntrons.hs new file mode 100644 --- /dev/null +++ b/test/GtfIntrons.hs @@ -0,0 +1,41 @@ +{-# LANGUAGE OverloadedStrings #-} +module Main + where + +import Control.Monad +import qualified Data.ByteString.Char8 as BS +import Data.List +import Data.Maybe +import System.Environment +import System.IO + +import Bio.SeqLoc.GTF +import Bio.SeqLoc.LocRepr +import qualified Bio.SeqLoc.Location as Loc +import Bio.SeqLoc.OnSeq +import qualified Bio.SeqLoc.SpliceLocation as SpLoc +import Bio.SeqLoc.Transcript +import Bio.SeqLoc.TranscriptTable + +main :: IO () +main = getArgs >>= mainWithArgs + where mainWithArgs [gtf] = gtfIntrons gtf + mainWithArgs _ = do prog <- getProgName + hPutStrLn stderr . unwords $ [ "USAGE:", prog, "" ] + +gtfIntrons :: FilePath -> IO () +gtfIntrons = readGtfTranscripts >=> mapM_ BS.putStr . intronGtf + where intronGtf = transcriptsToGtf . concatMap transcriptIntrons + +transcriptsToGtf :: [Transcript] -> [BS.ByteString] +transcriptsToGtf = map (transcriptToGtf "exons-to-introns") + +transcriptIntrons :: Transcript -> [Transcript] +transcriptIntrons trx = zipWith intronTranscript [1..] . junctions $ sploc + where (OnSeq refname sploc) = location trx + intronTranscript idx jct = Transcript geneid trxid loc Nothing + where geneid = SeqName . flip BS.append "_introns" . unSeqName . geneId $ trx + trxid = SeqName . flip BS.append trxsuffix . unSeqName . trxId $ trx + trxsuffix = "_intron" `BS.append` (BS.pack . show $ idx) + loc = OnSeq refname (fromMaybe noLoc $ SpLoc.fromContigs [ intron jct ]) + noLoc = error $ "Unable to create singleton SpLoc from " ++ (BS.unpack . repr) jct diff --git a/test/TestBed.hs b/test/TestBed.hs new file mode 100644 --- /dev/null +++ b/test/TestBed.hs @@ -0,0 +1,32 @@ +{-# LANGUAGE OverloadedStrings #-} +module Main + where + +import Control.Monad +import qualified Data.ByteString.Char8 as BS +import Data.List +import System.IO + +import qualified Data.Iteratee as Iter + +import Bio.SeqLoc.Bed +import Bio.SeqLoc.LocRepr +import Bio.SeqLoc.OnSeq +import Bio.SeqLoc.Transcript +import Bio.SeqLoc.TranscriptTable + +main :: IO () +main = do withFile "test/bed-out.txt" WriteMode $ \hout -> + let bedIter = bedTranscriptEnum $ Iter.mapM_ (BS.hPutStrLn hout . unparseLine) + in Iter.fileDriver bedIter "/data/genomes/Homo_sapiens/hg19_knownGene.bed" + withFile "test/bed-copy.bed" WriteMode $ \hout -> + let bedIter = bedTranscriptEnum $ Iter.mapM_ (BS.hPutStrLn hout . transcriptToBedStd) + in Iter.fileDriver bedIter "/data/genomes/Homo_sapiens/hg19_knownGene.bed" + + -- trx' <- readTable "test/out.txt" + -- let trx10 = take 10 trx' + -- BS.writeFile "test/out10.gtf" . BS.concat . map (transcriptToGtf "TestGtf") $ trx10 + -- trx10' <- readGtfTranscripts "test/out10.gtf" + -- print $ (sort . map location $ trx10) == (sort . map location $ trx10') + -- return () + diff --git a/test/TestGtf.hs b/test/TestGtf.hs new file mode 100644 --- /dev/null +++ b/test/TestGtf.hs @@ -0,0 +1,25 @@ +{-# LANGUAGE OverloadedStrings #-} +module Main + where + +import Control.Monad +import qualified Data.ByteString.Char8 as BS +import Data.List +import System.IO + +import Bio.SeqLoc.GTF +import Bio.SeqLoc.LocRepr +import Bio.SeqLoc.OnSeq +import Bio.SeqLoc.Transcript +import Bio.SeqLoc.TranscriptTable + +main :: IO () +main = do trx <- readGtfTranscripts "/data/genomes/Homo_sapiens/hg19_knownGene.gtf" + writeTable "test/gtf-out.txt" trx + trx' <- readTable "test/gtf-out.txt" + let trx10 = take 10 trx' + BS.writeFile "test/gtf-out10.gtf" . BS.concat . map (transcriptToGtf "TestGtf") $ trx10 + trx10' <- readGtfTranscripts "test/gtf-out10.gtf" + print $ (sort . map location $ trx10) == (sort . map location $ trx10') + return () +