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.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 diff --git a/src/BedSubregion.hs b/src/BedSubregion.hs new file mode 100644 --- /dev/null +++ b/src/BedSubregion.hs @@ -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 diff --git a/src/Bio/SeqLoc/Bed.hs b/src/Bio/SeqLoc/Bed.hs --- a/src/Bio/SeqLoc/Bed.hs +++ b/src/Bio/SeqLoc/Bed.hs @@ -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. diff --git a/src/GenomeToTrx.hs b/src/GenomeToTrx.hs new file mode 100644 --- /dev/null +++ b/src/GenomeToTrx.hs @@ -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 + + + + + diff --git a/test/GtfIntrons.hs b/test/GtfIntrons.hs --- a/test/GtfIntrons.hs +++ b/test/GtfIntrons.hs @@ -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, "" ] - -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)