seqloc-datafiles 0.2.2.2 → 0.4
raw patch · 5 files changed
+424/−25 lines, 5 filesdep +cmdthelinedep +conduitdep +conduit-extradep ~seqlocnew-component:exe:bed-subregionnew-component:exe:genome-to-trx
Dependencies added: cmdtheline, conduit, conduit-extra, filepath, lifted-base, pretty, resourcet, transformers-base, vector
Dependency ranges changed: seqloc
Files
- seqloc-datafiles.cabal +32/−11
- src/BedSubregion.hs +189/−0
- src/Bio/SeqLoc/Bed.hs +16/−1
- src/GenomeToTrx.hs +145/−0
- test/GtfIntrons.hs +42/−13
seqloc-datafiles.cabal view
@@ -1,5 +1,5 @@ Name: seqloc-datafiles-Version: 0.2.2.2+Version: 0.4 Cabal-Version: >= 1.6 Synopsis: Read and write BED and GTF format genome annotations Description: Read and write BED and GTF format genome annotations@@ -25,6 +25,7 @@ Other-modules: Bio.SeqLoc.ZeptoUtils Build-depends: base >= 4.2 && < 5, bytestring, attoparsec >= 0.8.5, hashable, unordered-containers,+ lifted-base >= 0.2.3, transformers-base >= 0.4.3, iteratee >= 0.8.1, seqloc >= 0.3.1, biocore >= 0.2 Hs-Source-Dirs: src Ghc-options: -Wall@@ -49,6 +50,28 @@ Ghc-options: -Wall -rtsopts C-Sources: src/rtsopts.c +Executable gtf-introns+ Main-is: GtfIntrons.hs+ Other-modules: Bio.SeqLoc.GTF, Bio.SeqLoc.TranscriptTable, Bio.SeqLoc.ZeptoUtils+ Build-depends: base >= 4.2 && < 5, bytestring,+ attoparsec >= 0.8.5, hashable, unordered-containers,+ iteratee >= 0.8.1, seqloc >= 0.3.1, biocore >= 0.2,+ QuickCheck, random, cmdtheline+ Hs-Source-Dirs: src, test+ Ghc-options: -Wall -rtsopts+ C-Sources: src/rtsopts.c++Executable bed-subregion+ Main-is: BedSubregion.hs+ Other-modules: Bio.SeqLoc.Bed, Bio.SeqLoc.ZeptoUtils+ Build-depends: base >= 4.2 && < 5, bytestring,+ attoparsec >= 0.8.5, hashable, unordered-containers,+ iteratee >= 0.8.1, seqloc >= 0.3.1, biocore >= 0.2, transformers, monads-tf,+ pretty+ Hs-Source-Dirs: src+ Ghc-options: -Wall -rtsopts+ C-Sources: src/rtsopts.c+ Executable test-gtf if !flag(Tests) Buildable: False@@ -74,15 +97,13 @@ 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+Executable genome-to-trx+ Main-is: GenomeToTrx.hs+ Buildable: False+ Other-modules: Bio.SeqLoc.Bed Build-depends: base >= 4.2 && < 5, bytestring, attoparsec >= 0.8.5, hashable, unordered-containers,- iteratee >= 0.8.1, seqloc >= 0.3.1, biocore >= 0.2,- QuickCheck, random- Hs-Source-Dirs: src, test- Ghc-options: -Wall -rtsopts- C-Sources: src/rtsopts.c+ iteratee >= 0.8.1, seqloc >= 0.6, biocore >= 0.2,+ conduit, conduit-extra, vector, filepath, resourcet+ Hs-Source-Dirs: src+ Ghc-options: -Wall
+ src/BedSubregion.hs view
@@ -0,0 +1,189 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE FlexibleContexts #-}+module Main+ where ++import Control.Applicative+import qualified Control.Exception.Lifted as E+import Control.Monad+import Control.Monad.Base+import Control.Monad.IO.Class+import qualified Control.Monad.Trans.Resource as R+import qualified Data.ByteString.Char8 as BS+import qualified Data.Conduit as C+import qualified Data.Conduit.Binary as CB+import qualified Data.Conduit.List as C+import Data.Maybe++import qualified Bio.SeqLoc.Bed as Bed+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 Loc+import Bio.SeqLoc.Strand+import Bio.SeqLoc.Transcript++import System.Console.CmdTheLine+import System.FilePath+import System.IO+import qualified Text.PrettyPrint as PP++main :: IO ()+main = run ( bedsub, info )+ where info = defTI { termName = "bed-subregion"+ , version = "2014-11-23"+ , termDoc = "Process BED file to extract specific sub-regions"+ }+ bedsub = bedSubregions <$> argConf++bedSubregions :: Conf -> IO ()+bedSubregions conf = R.runResourceT $+ CB.sourceFile (cInput conf) C.$=+ Bed.bedConduit C.$$+ handleTranscripts conf++handleTranscripts :: (MonadIO m, R.MonadResource m) => Conf -> C.Sink Transcript m ()+handleTranscripts conf = C.bracketP (openFile (cOutFile conf) WriteMode) hClose loop+ where loop hout = C.head >>= maybe (return ()) (\t -> handleOne t >> loop hout)+ where handleOne t = maybe (return ()) writeSubregion $ regionSpliceLoc (cRegionSpec conf) t+ where writeSubregion sl = let t' = t { location = (location t) { unOnSeq = sl }, cds = Nothing }+ in liftIO $ BS.hPutStrLn hout $ Bed.transcriptToBedStd t'+ +data TranscriptRegion = WholeTrx | Utr5 | Cds | Utr3 deriving (Show, Read, Eq, Ord, Bounded, Enum)++trxRegion :: TranscriptRegion -> Transcript -> Maybe Loc.ContigLoc+trxRegion WholeTrx trx = let sploc = unOnSeq . location $ trx+ in Just $! Loc.fromPosLen (Pos.Pos 0 Plus) (Loc.length sploc)+trxRegion Utr5 trx = utr5 trx+trxRegion Utr3 trx = utr3 trx+trxRegion Cds trx = cds trx+ +data TooLong = TooLongExtend | TooLongTruncate | TooLongDiscard deriving (Show, Read, Eq, Ord, Bounded, Enum)++handleEnds :: TooLong -> Loc.ContigLoc -> Loc.ContigLoc -> Maybe Loc.ContigLoc+handleEnds TooLongExtend _base cloc = Just cloc+handleEnds TooLongTruncate base cloc = let start = max (Pos.offset . Loc.startPos $ base) (Pos.offset . Loc.startPos $ cloc)+ end = min (Pos.offset . Loc.endPos $ base) (Pos.offset . Loc.endPos $ cloc)+ in if start < end+ then Just $! Loc.fromStartEnd start end+ else Nothing+handleEnds TooLongDiscard base cloc = if ((Pos.offset . Loc.startPos $ base) <= (Pos.offset . Loc.startPos $ cloc)+ && (Pos.offset . Loc.endPos $ base) >= (Pos.offset . Loc.endPos $ cloc))+ then Just cloc+ else Nothing++clocOutofExtended :: (Loc.Location l, LocRepr l) => Loc.ContigLoc -> l -> l+clocOutofExtended cloc l0 = fromMaybe badExtension $ Loc.clocOutof (Loc.slide startxtn cloc) lext+ where startxtn = negate . Pos.offset . Loc.startPos $ cloc+ endxtn = (Pos.offset . Loc.endPos $ cloc) + 1 - (Loc.length l0 - 1)+ lext = Loc.extend (startxtn, endxtn) l0+ badExtension = error . BS.unpack $ BS.unwords [ "Bad extension: ", repr cloc, " in ", repr l0, " to ", repr lext ]++data RegionSpec = RegionSpec { rsRegion :: !TranscriptRegion+ , rsStartOffset :: !(Maybe Pos.Offset)+ , rsLength :: !(Maybe Pos.Offset)+ , rsEndOffset :: !(Maybe Pos.Offset)+ , rsTooLong :: !TooLong+ }+ deriving (Show)++validRegionSpec :: RegionSpec -> Bool+validRegionSpec (RegionSpec _rgn (Just _startoff) (Just _len) Nothing _toolong) = True+validRegionSpec (RegionSpec _rgn (Just _startoff) Nothing (Just _endoff) _toolong) = True+validRegionSpec (RegionSpec _rgn Nothing (Just _len) (Just _endoff) _toolong) = True+validRegionSpec _ = False++regionSpliceLoc :: RegionSpec -> Transcript -> Maybe Loc.SpliceLoc+regionSpliceLoc (RegionSpec rgn (Just startoff) (Just len) Nothing toolong) trx+ = do base <- trxRegion rgn trx+ cloc <- handleEnds toolong base $ Loc.fromPosLen (Pos.slide (Loc.startPos base) startoff) len+ return $! clocOutofExtended cloc (unOnSeq . location $ trx)+regionSpliceLoc (RegionSpec rgn (Just startoff) Nothing (Just endoff) toolong) trx+ = do base <- trxRegion rgn trx+ let start = (Pos.offset . Loc.startPos $ base) + startoff+ end = (Pos.offset . Loc.endPos $ base) + endoff+ if start <= end+ then do cloc <- handleEnds toolong base $ Loc.fromStartEnd start end+ return $! clocOutofExtended cloc (unOnSeq . location $ trx)+ else Nothing+regionSpliceLoc (RegionSpec rgn Nothing (Just len) (Just endoff) toolong) trx+ = do base <- trxRegion rgn trx+ let end = (Pos.offset . Loc.endPos $ base) + endoff+ start = 1 + end - len+ if start <= end+ then do cloc <- handleEnds toolong base $ Loc.fromStartEnd start end+ return $! clocOutofExtended cloc (unOnSeq . location $ trx)+ else Nothing+regionSpliceLoc rs _ = error $ "Invalid region selection " ++ show rs++data Conf = Conf { cInput :: !FilePath+ , cOutput :: !(Maybe FilePath)+ , cRegionSpec :: !RegionSpec+ }++cOutFile :: Conf -> FilePath+cOutFile conf = fromMaybe defaultOutput . cOutput $ conf+ where defaultOutput = (dropExtension . cInput $ conf) ++ "_subregion" ++ (takeExtension . cInput $ conf)++argConf :: Term Conf+argConf = Conf <$> argInput <*> argOutput <*> regionspec++argInput :: Term FilePath+argInput = required $ opt Nothing $ ( optInfo [ "i" ])+ { optName = "INPUT.BED", optDoc = "BED input" }++argOutput :: Term (Maybe FilePath)+argOutput = value $ opt Nothing $ ( optInfo [ "o" ])+ { optName = "OUTPUT.BED", optDoc = "BED output" }++argRegion :: Term TranscriptRegion+argRegion = ret . fmap validate . value $+ vFlag Nothing $+ [ (Just WholeTrx, (optInfo [ "whole-trx" ]) { optDoc = "Region relative to whole transcript" })+ , (Just Utr5, (optInfo [ "utr5" ]) { optDoc = "Region relative to 5' UTR" })+ , (Just Utr3, (optInfo [ "utr3" ]) { optDoc = "Region relative to 3' UTR" })+ , (Just Cds, (optInfo [ "cds" ]) { optDoc = "Region relative to CDS" })+ ]+ where validate :: Maybe TranscriptRegion -> Err TranscriptRegion+ validate = maybe noRegion return+ noRegion = msgFail . PP.text $+ "Specify a reference transcription region (whole transcript, CDS, etc.)"++argTooLong :: Term TooLong+argTooLong = ret . fmap validate . value $+ vFlag Nothing $+ [ (Just TooLongExtend, (optInfo [ "extend" ]) { optDoc = "Extend beyond reference region" })+ , (Just TooLongTruncate, (optInfo [ "truncate" ]) { optDoc = "Truncate to lie within reference region " })+ , (Just TooLongDiscard, (optInfo [ "discard" ]) { optDoc = "Discard when lying outside reference region" })+ ]+ where validate :: Maybe TooLong -> Err TooLong+ validate = maybe noarg return+ noarg = msgFail . PP.text $+ "Specify how subregions extending outside the reference region should be handled (truncation etc.)"++argStartOffset :: Term (Maybe Int)+argStartOffset = value $ opt Nothing $ ( optInfo [ "start-off" ]) { optName = "DELTA-START", optDoc = "Offset of start position, positive is more 3'" }++argEndOffset :: Term (Maybe Int)+argEndOffset = value $ opt Nothing $ ( optInfo [ "end-off" ]) { optName = "DELTA-END", optDoc = "Offset of end position, positive is more 3'" }++argLength :: Term (Maybe Int)+argLength = value $ opt Nothing $ ( optInfo [ "length" ]) { optName = "LENGTH", optDoc = "Length of the sub-region" }++regionspec :: Term RegionSpec+regionspec = ret . fmap validate $+ RegionSpec <$>+ argRegion <*>+ (liftM fromIntegral <$> argStartOffset) <*>+ (liftM fromIntegral <$> argLength) <*>+ (liftM fromIntegral <$> argEndOffset) <*>+ argTooLong+ where validate :: RegionSpec -> Err RegionSpec+ validate rs = if validRegionSpec rs+ then return rs+ else msgFail . PP.text $+ "Specify exactly two of start offset, length, and end offset"++throwerr :: (MonadBase IO m) => String -> m a+throwerr = E.ioError . userError
src/Bio/SeqLoc/Bed.hs view
@@ -1,21 +1,27 @@-{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE OverloadedStrings, FlexibleContexts #-} {-| Utilities for reading and writing BED format gene annotations -} module Bio.SeqLoc.Bed ( readBedTranscripts , bedZP, bedTranscriptEnum+ , bedConduit , transcriptToBed, transcriptToBedStd ) where import Control.Applicative+import qualified Control.Exception.Lifted as E import Control.Monad+import Control.Monad.Base import qualified Data.ByteString.Char8 as BS import Data.List import Data.Maybe import Data.Ord import qualified Data.Attoparsec.Zepto as ZP+import qualified Data.Conduit as C+import qualified Data.Conduit.Binary as CB+import qualified Data.Conduit.List as C import qualified Data.Iteratee as Iter import qualified Data.Iteratee.Char as IterChar @@ -77,6 +83,15 @@ 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 ++-- | Conduit from a 'BS.ByteString' source such as a BED file to a+-- source of 'Transcript' annotations from the file.+bedConduit :: (Monad m, MonadBase IO m) => C.Conduit BS.ByteString m Transcript+bedConduit = CB.lines C.$= loop+ where loop = C.head >>= maybe (return ())+ (\l -> case ZP.parse bedZP l of+ Left err -> E.ioError . userError $ err ++ "\n in BED line\n" ++ show l+ Right res -> C.yield res >> loop) -- | Minimalistic 'ZP.Parser'-style parser for a BED format line, not -- including the trailing newline.
+ src/GenomeToTrx.hs view
@@ -0,0 +1,145 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Main+ where++import Control.Applicative+import Control.Monad+import qualified Data.ByteString.Char8 as BS+import Data.Char+import qualified Data.HashMap.Strict as HM+import Data.Maybe+import Numeric+import System.Exit+import System.FilePath+import System.IO++import qualified Data.Attoparsec.Zepto as ZP+import qualified Data.Vector.Unboxed as U+import System.Console.CmdTheLine++import qualified Control.Monad.Trans.Resource as R+import qualified Data.Conduit as C+import qualified Data.Conduit.Binary as CB+import qualified Data.Conduit.List as C++import qualified Bio.SeqLoc.Bed as Bed+import qualified Bio.SeqLoc.LocMap as LM+import qualified Bio.SeqLoc.Location as Loc+import Bio.SeqLoc.OnSeq+import qualified Bio.SeqLoc.Position as Pos+import Bio.SeqLoc.Strand+import Bio.SeqLoc.Transcript+import Bio.SeqLoc.ZeptoUtils++main :: IO ()+main = run ( g2t, info )+ where info = defTI { termName = "genome-to-trx"+ , version = "0.0"+ , termDoc = "Convert genome coordinates to transcriptome coordinates"+ }+ g2t = genomeToTrx <$> trxBed <*> coordInput <*> outputFile <*> argConf++genomeToTrx :: String -> String -> Maybe String -> Conf -> IO ()+genomeToTrx trx input moutput conf+ = do trxs <- liftM (LM.transcriptSeqLocMap 100000) $ Bed.readBedTranscripts trx+ R.runResourceT $ CB.sourceFile input C.$= remapCoordLines trxs conf C.$$ CB.sinkFile output+ where output = fromMaybe defaultOutput moutput+ defaultOutput = let (base, ext) = splitExtension input+ in (base ++ "_trx") <.> ext++data InputType = InputBed | InputBedPE deriving (Read, Show, Ord, Eq)+data Strandedness = FwdOnly | RevOnly | Both deriving (Read, Show, Ord, Eq)+data Coding = TrxRelative | CdsRelative | CdsOnly deriving (Read, Show, Ord, Eq)++data Conf = Conf { cInputType :: !InputType,+ cStrandedness :: !Strandedness,+ cCodingRelative :: !Coding,+ cReportNoHit :: !Bool+ } deriving (Read, Show)++trxBed :: Term String+trxBed = required $ opt Nothing $ (optInfo [ "t", "transcripts" ])+ { optName = "TRANSCRIPTS.BED", optDoc = "Transcript BED file" }++coordInput :: Term String+coordInput = required $ opt Nothing $ (optInfo [ "i", "input" ])+ { optName = "INPUT.TXT", optDoc = "Input coordinates" }++outputFile :: Term (Maybe String)+outputFile = value $ opt Nothing $ (optInfo [ "o", "output" ])+ { optName = "OUTPUT.TXT", optDoc = "Output coordinates" }++argConf :: Term Conf+argConf = Conf <$>+ confInputType <*>+ confStrandedness <*>+ confCodingRelative <*>+ confReportNoHit+ where confInputType :: Term InputType+ confInputType = value $ vFlag InputBed [(InputBedPE, (optInfo [ "bedpe" ]) { optName = "Bed PE (bedtools) format", optDoc = "Bed PE (bedtools) format" })]+ confStrandedness :: Term Strandedness+ confStrandedness = value $ vFlag FwdOnly+ [( FwdOnly, (optInfo [ "fwd" ]) { optName = "Forward feature strand only", optDoc = "Forward feature strand only" }),+ ( RevOnly, (optInfo [ "rev" ]) { optName = "Reverse feature strand only", optDoc = "Reverse feature strand only" }),+ ( Both, (optInfo [ "both" ]) { optName = "Both feature strands", optDoc = "Both feature strands" })]+ confCodingRelative :: Term Bool+ confCodingRelative = value $ vFlag Transcript + flag $ (optInfo [ "c", "cds" ]) { optName = "Coords relative to CDS start", optDoc = "Coords relative to CDS start" }+ confReportNoHit :: Term Bool+ confReportNoHit = value $ flag $ (optInfo [ "n", "no-hit" ]) { optName = "Report no-hit lines", optDoc = "Report no-hit lines" }++remapCoordLines :: (Monad m) => LM.SeqLocMap Transcript -> Conf -> C.Conduit BS.ByteString m BS.ByteString+remapCoordLines trxs conf = CB.lines C.=$= C.concatMapM remapCoordLine C.=$= C.map (flip BS.append "\n")+ where remapCoordLine :: (Monad m) => BS.ByteString -> m [BS.ByteString]+ remapCoordLine = case cInputType conf of+ InputBed -> remapBedLine trxs conf+ etc -> fail $ "Input type " ++ show etc ++ " not implemented"++remapBedLine :: (Monad m) => LM.SeqLocMap Transcript -> Conf -> BS.ByteString -> m [BS.ByteString]+remapBedLine trxs conf l = case BS.split '\t' l of+ (chr:startBS:endBS:name:score:strandBS:rest)+ -> do start <- either (const . fail $ "Bad start in " ++ show l) return $ ZP.parse decimal startBS+ end <- either (const . fail $ "Bad end in " ++ show l) return $ ZP.parse decimal endBS+ strand <- case strandBS of+ "+" -> return Plus+ "-" -> return Minus+ _ -> fail $ "Bad strand " ++ show strandBS ++ " in " ++ show l+ let gloc = OnSeq (toSeqLabel chr) (Loc.fromBoundsStrand start (end - 1) strand)+ let tlocs = remapLoc trxs conf gloc+ bedLine tsloc = BS.intercalate "\t" [ unSeqLabel . onSeqLabel $ tsloc+ , BS.pack . show . Pos.unOff . fst . Loc.bounds . unOnSeq $ tsloc+ , BS.pack . show . (+ 1) . Pos.unOff . snd . Loc.bounds . unOnSeq $ tsloc+ , name, score+ , if ((Loc.strand . unOnSeq) tsloc == Plus) then "+" else "-" ]+ if null tlocs && cReportNoHit conf+ then return $! [BS.intercalate "\t" $ [ "N/A", ".", ".", name, score, "." ]]+ else return $! map bedLine tlocs++remapLoc :: LM.SeqLocMap Transcript -> Conf -> ContigSeqLoc -> [ContigSeqLoc]+remapLoc trxs conf l = mapMaybe locInto cands+ where candList = LM.querySeqLoc l trxs+ cands = HM.elems . HM.fromList . map (\t -> (trxId t, t)) $ candList+ locInto t = let tsloc = location t+ in Loc.clocInto (unOnSeq l) (unOnSeq tsloc) >>= \tloc ->+ if (Loc.strand tloc == Plus && cStrandedness conf == RevOnly) ||+ (Loc.strand tloc == Minus && cStrandedness conf == FwdOnly)+ then Nothing+ else Just (OnSeq (trxId t) tloc)++locInto :: Transcript -> Conf -> ContigSeqLoc -> ContigSeqLoc+locinto t conf l = Loc.clocInto (unOnSeq l) (unOnSeq tsloc) >>=+ strandCheck >>= \tloc ->+ cdsAdjust (OnSeq (trxId t) tloc)+ where tsloc = location t+ strandCheck tloc = case cStrandedness conf of+ FwdOnly -> if (Loc.strand tloc == Plus) then Just tloc else Nothing+ RevOnly -> if (Loc.strand tloc == Minus) then Just tloc else Nothing+ Both -> Just tloc+ cdsAdjust = case cCodingRelative conf of+ + + + +
test/GtfIntrons.hs view
@@ -2,6 +2,7 @@ module Main where +import Control.Applicative import Control.Monad import qualified Data.ByteString.Char8 as BS import Data.List@@ -9,6 +10,8 @@ import System.Environment import System.IO +import System.Console.CmdTheLine+ import Bio.SeqLoc.GTF import Bio.SeqLoc.LocRepr import qualified Bio.SeqLoc.Location as Loc@@ -18,24 +21,50 @@ 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+main = run ( gtfIntrons, info )+ where info = defTI { termName = "gtf-introns"+ , version = "0.0"+ , termDoc = "Generate GTF file of introns from GTF file of transcripts"+ }+ gtfIntrons = introns <$> naming <*> gtfin+ introns nam = readGtfTranscripts >=> mapM_ BS.putStr . intronGtf+ where intronGtf = transcriptsToGtf . concatMap (transcriptIntrons nam) transcriptsToGtf :: [Transcript] -> [BS.ByteString] transcriptsToGtf = map (transcriptToGtf "exons-to-introns") -transcriptIntrons :: Transcript -> [Transcript]-transcriptIntrons trx = zipWith intronTranscript [1..] . junctions $ sploc+transcriptIntrons :: Naming -> Transcript -> [Transcript]+transcriptIntrons nam trx = zipWith intronTranscript [1..] . junctions $ sploc where (OnSeq refname sploc) = location trx intronTranscript idx jct = Transcript geneid trxid loc Nothing- where geneid = toSeqLabel . flip BS.append "_introns" . unSeqLabel . geneId $ trx- trxid = toSeqLabel . flip BS.append trxsuffix . unSeqLabel . trxId $ trx- trxsuffix = "_intron" `BS.append` (BS.pack . show $ idx)+ where geneid = intronGeneId nam idx trx -- toSeqLabel . flip BS.append "_introns" . unSeqLabel . geneId $ trx+ trxid = intronTrxId nam idx trx -- toSeqLabel . flip BS.append trxsuffix . unSeqLabel . trxId $ trx loc = OnSeq refname (fromMaybe noLoc $ SpLoc.fromContigs [ intron jct ]) noLoc = error $ "Unable to create singleton SpLoc from " ++ (BS.unpack . repr) jct+ -- trxsuffix = "_intron" `BS.append` (BS.pack . show $ idx) + +naming :: Term Naming+naming = value $ vFlag SameGene + [ ( DifferentGenes, (optInfo [ "different-genes" ]) { optDoc = "Different GTF genes for each intron" } )+ , ( SameGene, (optInfo [ "same-gene" ]) { optDoc = "Different GTF transcripts (but the same GTF gene) for each intron" })+ , ( SameTranscript, (optInfo [ "same-transcript" ]) { optDoc = "Same GTF transcript for each intron" } )+ ]++gtfin :: Term String+gtfin = required $ pos 0 Nothing $ posInfo { posName = "GTF", posDoc = "Input GTF file" }++data Naming = DifferentGenes+ | SameGene+ | SameTranscript+ +intronGeneId :: Naming -> Int -> Transcript -> SeqLabel+intronGeneId naming idx trx = toSeqLabel . flip BS.append genesuffix . unSeqLabel . geneId $ trx+ where genesuffix = case naming of+ DifferentGenes -> "_intron" `BS.append` (BS.pack . show $ idx)+ _ -> "_introns"+ +intronTrxId :: Naming -> Int -> Transcript -> SeqLabel+intronTrxId naming idx trx = toSeqLabel . flip BS.append trxsuffix . unSeqLabel . geneId $ trx+ where trxsuffix = case naming of+ SameTranscript -> "_introns"+ _ -> "_intron" `BS.append` (BS.pack . show $ idx)