seqloc-datafiles (empty) → 0.1
raw patch · 12 files changed
+657/−0 lines, 12 filesdep +QuickCheckdep +attoparsecdep +basesetup-changed
Dependencies added: QuickCheck, attoparsec, base, bytestring, hashable, haskell98, iteratee, monads-tf, random, seqloc, transformers, unordered-containers
Files
- LICENSE +21/−0
- Setup.hs +2/−0
- seqloc-datafiles.cabal +73/−0
- src/Bio/SeqLoc/Bed.hs +144/−0
- src/Bio/SeqLoc/GTF.hs +187/−0
- src/Bio/SeqLoc/TranscriptTable.hs +48/−0
- src/Bio/SeqLoc/ZeptoUtils.hs +32/−0
- src/GtfToBed.hs +51/−0
- src/rtsopts.c +1/−0
- test/GtfIntrons.hs +41/−0
- test/TestBed.hs +32/−0
- test/TestGtf.hs +25/−0
+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ seqloc-datafiles.cabal view
@@ -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
+ src/Bio/SeqLoc/Bed.hs view
@@ -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)
+ src/Bio/SeqLoc/GTF.hs view
@@ -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 "\";"+ ++
+ src/Bio/SeqLoc/TranscriptTable.hs view
@@ -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"
+ src/Bio/SeqLoc/ZeptoUtils.hs view
@@ -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)
+ src/GtfToBed.hs view
@@ -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
+ src/rtsopts.c view
@@ -0,0 +1,1 @@+char *ghc_rts_opts = "-A128m";
+ test/GtfIntrons.hs view
@@ -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, "<GTF>" ]+ +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
+ test/TestBed.hs view
@@ -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 ()+
+ test/TestGtf.hs view
@@ -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 ()+