diff --git a/seqloc-datafiles.cabal b/seqloc-datafiles.cabal --- a/seqloc-datafiles.cabal +++ b/seqloc-datafiles.cabal @@ -1,6 +1,6 @@ Name: seqloc-datafiles -Version: 0.4 -Cabal-Version: >= 1.6 +Version: 0.4.2 +Cabal-Version: >= 1.10 Synopsis: Read and write BED and GTF format genome annotations Description: Read and write BED and GTF format genome annotations License: MIT @@ -16,39 +16,41 @@ location: http://github.com/ingolia/seqloc subdir: seqloc-datafiles -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.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 + iteratee >= 0.8.1, seqloc >= 0.3.1, biocore >= 0.2, + conduit, conduit-extra, resourcet Hs-Source-Dirs: src Ghc-options: -Wall + Default-Language: Haskell2010 Executable gtf-to-bed Main-is: GtfToBed.hs Other-modules: Bio.SeqLoc.GTF, 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 + iteratee >= 0.8.1, seqloc >= 0.3.1, biocore >= 0.2, transformers, monads-tf, + lifted-base, transformers-base, conduit, conduit-extra, resourcet Hs-Source-Dirs: src Ghc-options: -Wall -rtsopts C-Sources: src/rtsopts.c + Default-Language: Haskell2010 Executable bed-to-gtf Main-is: BedToGtf.hs Other-modules: Bio.SeqLoc.GTF, 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 + iteratee >= 0.8.1, seqloc >= 0.3.1, biocore >= 0.2, transformers, monads-tf, + conduit, conduit-extra, lifted-base, transformers-base, resourcet Hs-Source-Dirs: src Ghc-options: -Wall -rtsopts C-Sources: src/rtsopts.c + Default-Language: Haskell2010 Executable gtf-introns Main-is: GtfIntrons.hs @@ -56,10 +58,12 @@ 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 + QuickCheck, random, cmdtheline, + conduit, conduit-extra, resourcet Hs-Source-Dirs: src, test Ghc-options: -Wall -rtsopts C-Sources: src/rtsopts.c + Default-Language: Haskell2010 Executable bed-subregion Main-is: BedSubregion.hs @@ -67,14 +71,27 @@ 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 + pretty, conduit, conduit-extra, resourcet, cmdtheline, filepath, lifted-base, + transformers-base, resourcet Hs-Source-Dirs: src Ghc-options: -Wall -rtsopts C-Sources: src/rtsopts.c + Default-Language: Haskell2010 -Executable test-gtf - if !flag(Tests) - Buildable: False +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.6, biocore >= 0.2, + conduit, conduit-extra, vector, filepath, resourcet + Hs-Source-Dirs: src + Ghc-options: -Wall + Default-Language: Haskell2010 + +Test-suite test-gtf + Type: exitcode-stdio-1.0 Main-is: TestGtf.hs Other-modules: Bio.SeqLoc.GTF, Bio.SeqLoc.TranscriptTable, Bio.SeqLoc.ZeptoUtils Build-depends: base >= 4.2 && < 5, bytestring, @@ -84,26 +101,28 @@ Hs-Source-Dirs: src, test Ghc-options: -Wall -rtsopts C-Sources: src/rtsopts.c + Default-Language: Haskell2010 -Executable test-bed - if !flag(Tests) - Buildable: False +Test-Suite test-bed + Type: exitcode-stdio-1.0 Main-is: TestBed.hs Other-modules: Bio.SeqLoc.Bed, 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 + QuickCheck, random, lifted-base, transformers-base, + conduit, conduit-extra Hs-Source-Dirs: src, test Ghc-options: -Wall + Default-Language: Haskell2010 -Executable genome-to-trx - Main-is: GenomeToTrx.hs - Buildable: False +Test-Suite test-subregions + Type: exitcode-stdio-1.0 + Main-is: TestSubregions.hs 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.6, biocore >= 0.2, - conduit, conduit-extra, vector, filepath, resourcet - Hs-Source-Dirs: src - Ghc-options: -Wall + Build-depends: base >= 4.2 && < 5, bytestring, attoparsec >= 0.8.5, + QuickCheck, random, biocore >= 0.2, unordered-containers, vector, hashable, + seqloc >= 0.5, lifted-base, transformers-base, conduit, conduit-extra, + iteratee, directory, process, transformers + Hs-Source-Dirs: src, test + Default-Language: Haskell2010 diff --git a/src/BedSubregion.hs b/src/BedSubregion.hs --- a/src/BedSubregion.hs +++ b/src/BedSubregion.hs @@ -13,6 +13,7 @@ import qualified Data.Conduit as C import qualified Data.Conduit.Binary as CB import qualified Data.Conduit.List as C +import Data.List import Data.Maybe import qualified Bio.SeqLoc.Bed as Bed @@ -47,11 +48,24 @@ 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' + where writeSubregion sl = liftIO . BS.hPutStrLn hout . Bed.transcriptToBedStd $ subtranscript sl + subtranscript sl = Transcript { geneId = gene', trxId = trx', location = loc', cds = Nothing } + where loc' = (location t) { unOnSeq = sl } + gene' = toSeqLabel . flip BS.append suffix . unSeqLabel . geneId $ t + trx' = toSeqLabel . flip BS.append suffix . unSeqLabel . trxId $ t + suffix = case cRenameFeatures conf of + NoRename -> BS.empty + AutoRename -> BS.pack . ('_' :) . regionSpecName . cRegionSpec $ conf + Rename sfx -> BS.pack sfx data TranscriptRegion = WholeTrx | Utr5 | Cds | Utr3 deriving (Show, Read, Eq, Ord, Bounded, Enum) +regionName :: TranscriptRegion -> String +regionName WholeTrx = "trx" +regionName Utr5 = "utr5" +regionName Cds = "cds" +regionName Utr3 = "utr3" + 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) @@ -61,12 +75,17 @@ data TooLong = TooLongExtend | TooLongTruncate | TooLongDiscard deriving (Show, Read, Eq, Ord, Bounded, Enum) +tooLongName :: TooLong -> String +tooLongName TooLongExtend = "ext" +tooLongName TooLongTruncate = "trunc" +tooLongName TooLongDiscard = "disc" + 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 + in if start <= end + then Just $! Loc.fromBoundsStrand start end Plus 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)) @@ -88,6 +107,21 @@ } deriving (Show) +regionSpecName :: RegionSpec -> String +regionSpecName (RegionSpec rgn (Just startoff) (Just len) Nothing toolong) + = intercalate "_" [ regionName rgn, "start" ++ showSigned startoff, "length" ++ (show . Pos.unOff $ len), tooLongName toolong ] +regionSpecName (RegionSpec rgn (Just startoff) Nothing (Just endoff) toolong) + = intercalate "_" [ regionName rgn, "start" ++ showSigned startoff, "end" ++ show endoff, tooLongName toolong ] +regionSpecName (RegionSpec rgn Nothing (Just len) (Just endoff) toolong) + = intercalate "_" [ regionName rgn, "length" ++ (show . Pos.unOff $ len), "end" ++ show endoff, tooLongName toolong ] +regionSpecName rs = error $ "regionSpecName: invalid region specification " ++ show rs + +showSigned :: Pos.Offset -> String +showSigned z = case show . Pos.unOff $ z of + [] -> [] + str@('-':_) -> str + str -> '+' : str + validRegionSpec :: RegionSpec -> Bool validRegionSpec (RegionSpec _rgn (Just _startoff) (Just _len) Nothing _toolong) = True validRegionSpec (RegionSpec _rgn (Just _startoff) Nothing (Just _endoff) _toolong) = True @@ -104,7 +138,7 @@ 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 + then do cloc <- handleEnds toolong base $ Loc.fromBoundsStrand start end Plus return $! clocOutofExtended cloc (unOnSeq . location $ trx) else Nothing regionSpliceLoc (RegionSpec rgn Nothing (Just len) (Just endoff) toolong) trx @@ -112,22 +146,29 @@ 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 + then do cloc <- handleEnds toolong base $ Loc.fromBoundsStrand start end Plus return $! clocOutofExtended cloc (unOnSeq . location $ trx) else Nothing regionSpliceLoc rs _ = error $ "Invalid region selection " ++ show rs +data Rename = NoRename | AutoRename | Rename { unRename :: !String } deriving (Show) + +instance ArgVal Rename where + converter = let (stringParser, stringPrinter) = converter + in (either Left (Right . Rename) . stringParser, stringPrinter . unRename) + data Conf = Conf { cInput :: !FilePath , cOutput :: !(Maybe FilePath) , cRegionSpec :: !RegionSpec - } + , cRenameFeatures :: !Rename + } deriving (Show) cOutFile :: Conf -> FilePath cOutFile conf = fromMaybe defaultOutput . cOutput $ conf - where defaultOutput = (dropExtension . cInput $ conf) ++ "_subregion" ++ (takeExtension . cInput $ conf) + where defaultOutput = (dropExtension . cInput $ conf) ++ "_" ++ regionSpecName (cRegionSpec conf) ++ (takeExtension . cInput $ conf) argConf :: Term Conf -argConf = Conf <$> argInput <*> argOutput <*> regionspec +argConf = Conf <$> argInput <*> argOutput <*> regionspec <*> argRename argInput :: Term FilePath argInput = required $ opt Nothing $ ( optInfo [ "i" ]) @@ -184,6 +225,10 @@ then return rs else msgFail . PP.text $ "Specify exactly two of start offset, length, and end offset" + +argRename :: Term Rename +argRename = value $ defaultOpt AutoRename NoRename $ ( optInfo [ "r", "rename-features" ]) + { optName = "SUFFIX", optDoc = "Rename features to describe subregion [defaults to systematic name based on region specifier]" } 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 @@ -4,7 +4,7 @@ module Bio.SeqLoc.Bed ( readBedTranscripts , bedZP, bedTranscriptEnum - , bedConduit + , bedConduit, unbedConduit , transcriptToBed, transcriptToBedStd ) where @@ -13,6 +13,7 @@ import qualified Control.Exception.Lifted as E import Control.Monad import Control.Monad.Base +import Control.Monad.Trans.Resource import qualified Data.ByteString.Char8 as BS import Data.List import Data.Maybe @@ -72,8 +73,9 @@ -- | Read all BED format annotations in a BED file readBedTranscripts :: FilePath -> IO [Transcript] -readBedTranscripts = Iter.fileDriver (bedTranscriptEnum Iter.stream2list) - +readBedTranscripts bedfile = runResourceT $ C.runConduit $ + CB.sourceFile bedfile C.$= bedConduit C.$= C.consume + -- | 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. @@ -93,6 +95,10 @@ Left err -> E.ioError . userError $ err ++ "\n in BED line\n" ++ show l Right res -> C.yield res >> loop) +unbedConduit :: (Monad m) => C.Conduit Transcript m BS.ByteString +unbedConduit = C.head >>= maybe (return ()) go + where go t = C.yield (transcriptToBedStd t `BS.append` "\n") >> unbedConduit + -- | Minimalistic 'ZP.Parser'-style parser for a BED format line, not -- including the trailing newline. bedZP :: ZP.Parser Transcript @@ -158,7 +164,13 @@ = maybe badCdsLoc return $ do relstart <- Loc.posInto (Pos.Pos thickStart Plus) loc relend <- Loc.posInto (Pos.Pos (thickEnd - 1) Plus) loc - return $! stranded (Loc.strand loc) $ Loc.fromStartEnd (Pos.offset relstart) (Pos.offset relend) + case Loc.strand loc of + Plus -> if relstart <= relend + then return $! Loc.fromBoundsStrand (Pos.offset relstart) (Pos.offset relend) Plus + else Nothing + Minus -> if relend <= relstart + then return $! Loc.fromBoundsStrand (Pos.offset relend) (Pos.offset relstart) Plus + else Nothing where badCdsLoc = fail $ "Bio.SeqLoc.Bed: bad cds in " ++ (BS.unpack . BS.unwords $ [ repr loc, repr thickStart, repr thickEnd ]) diff --git a/test/TestBed.hs b/test/TestBed.hs --- a/test/TestBed.hs +++ b/test/TestBed.hs @@ -16,12 +16,17 @@ 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" +main = do trxs <- readBedTranscripts "/mnt/ingolialab/ingolia/Genomes/Saccharomyces_cerevisiae/YeastGenome/sac_cer_yassour.bed" + hPutStrLn stderr $! "Got " ++ (show . length $ trxs) ++ " transcripts" + trxs <- readBedTranscripts "/mnt/ingolialab/ingolia/Genomes/Homo_sapiens/GRCh38/gencode.bed" + hPutStrLn stderr $! "Got " ++ (show . length $ trxs) ++ " transcripts" + +-- 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' diff --git a/test/TestSubregions.hs b/test/TestSubregions.hs new file mode 100644 --- /dev/null +++ b/test/TestSubregions.hs @@ -0,0 +1,214 @@ +{-# LANGUAGE BangPatterns, ExistentialQuantification #-} +module Main + where + +import Control.Applicative +import Control.Monad +import Control.Monad.IO.Class +import qualified Data.ByteString.Char8 as BS +import Data.List +import Data.Maybe +import System.Directory +import System.IO +import System.Process +import Test.QuickCheck +import Test.QuickCheck.Monadic + +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 SpLoc +import Bio.SeqLoc.Strand +import Bio.SeqLoc.Transcript + +bedSubregion = "./dist/build/bed-subregion/bed-subregion" + +main :: IO () +main = mapM_ runTest tests + +tests :: [Test] +tests = [ T "Write random bed" test_writeRandomBed + , T "Whole transcript" test_wholeTranscript + , T "CDS only" test_CDS + , T "CDS start & length" test_CDS_start_length + , T "Trx start & length" test_trx_start_length + , T "Trx start & end" test_trx_start_end + , T "Trx end & length" test_trx_end_length + , T "Truncate start" test_truncate_start + , T "Truncate end" test_truncate_end + ] + +test_writeRandomBed :: Property +test_writeRandomBed = monadicIO $ + forAllM (listOf1 genTranscript) $ \trxs -> + run . BS.writeFile "tmpf.bed" . BS.unlines . map Bed.transcriptToBedStd $ trxs + +test_wholeTranscript :: Property +test_wholeTranscript = test_subtranscript (Just . stripCDS) [ "--whole-trx", "--discard", "--start=0", "--end=0" ] + where stripCDS t = t { cds = Nothing } + +test_CDS :: Property +test_CDS = test_subtranscript subcds [ "--cds", "--discard", "--start=0", "--end=0" ] + where subcds t = let (OnSeq chr loc) = location t + in do subloc <- cds t >>= \cdsloc -> Loc.clocOutof cdsloc loc + return $! t { location = (OnSeq chr subloc), cds = Nothing } + +test_CDS_start_length :: Property +test_CDS_start_length = forAll genNonNegOffset $ \substart -> + forAll genPositiveOffset $ \sublen -> + let subcdslen t = let (OnSeq chr loc) = location t + subsubloc = Loc.fromPosLen (Pos.Pos substart Plus) sublen + in do subcdsloc <- cds t >>= Loc.clocOutof subsubloc + subloc <- Loc.clocOutof subcdsloc loc + return $! t { location = (OnSeq chr subloc), cds = Nothing } + in test_subtranscript subcdslen [ "--cds", "--discard", "--start=" ++ show (Pos.unOff substart), "--length=" ++ show (Pos.unOff sublen) ] + +test_trx_start_length :: Property +test_trx_start_length = forAll genNonNegOffset $ \substart -> + forAll genPositiveOffset $ \sublen -> + let subtrx t = let (OnSeq chr loc) = location t + subsubloc = Loc.fromBoundsStrand substart (substart + sublen - 1) Plus + in do subloc <- Loc.clocOutof subsubloc loc + return $! t { location = (OnSeq chr subloc), cds = Nothing } + in test_subtranscript subtrx [ "--whole-trx", "--discard", "--start=" ++ show (Pos.unOff substart), "--length=" ++ show (Pos.unOff sublen) ] + +test_trx_end_length :: Property +test_trx_end_length = forAll genNonNegOffset $ \subend -> + forAll genPositiveOffset $ \sublen -> + let subtrxlen t = let (OnSeq chr loc) = location t + trxlen = Loc.length loc + subsubloc = Loc.fromBoundsStrand (trxlen - (subend + sublen)) (trxlen - (subend + 1)) Plus + in do subloc <- Loc.clocOutof subsubloc loc + return $! t { location = (OnSeq chr subloc), cds = Nothing } + in test_subtranscript subtrxlen [ "--whole-trx", "--discard", "--end=" ++ show (negate $ Pos.unOff subend), "--length=" ++ show (Pos.unOff sublen) ] + +test_trx_start_end :: Property +test_trx_start_end = forAll genNonNegOffset $ \substart -> + forAll genNonNegOffset $ \subend -> + let subtrxlen t = let (OnSeq chr loc) = location t + trxlen = Loc.length loc + in do subsubloc <- if substart <= (trxlen - (subend + 1)) + then Just $! Loc.fromBoundsStrand substart (trxlen - (subend + 1)) Plus + else Nothing + subloc <- Loc.clocOutof subsubloc loc + return $! t { location = (OnSeq chr subloc), cds = Nothing } + in test_subtranscript subtrxlen [ "--whole-trx", "--discard", "--start=" ++ show (Pos.unOff substart), "--end=" ++ show (negate $ Pos.unOff subend) ] + +test_truncate_start :: Property +test_truncate_start = forAll genNonNegOffset $ \inlen -> + forAll genNonNegOffset $ \outlen -> + let subtrx t = let (OnSeq chr loc) = location t + trxlen = Loc.length loc + efflen = min inlen trxlen + in do subsubloc <- if efflen > 0 + then Just $! Loc.fromBoundsStrand 0 (efflen - 1) Plus + else Nothing + subloc <- Loc.clocOutof subsubloc loc + return $! t { location = (OnSeq chr subloc), cds = Nothing } + in test_subtranscript subtrx [ "--whole-trx", "--truncate", "--start=" ++ show (negate $ Pos.unOff outlen), "--length=" ++ show (Pos.unOff $ inlen + outlen) ] + +test_truncate_end :: Property +test_truncate_end = forAll genNonNegOffset $ \inlen -> + forAll genNonNegOffset $ \outlen -> + let subtrx t = let (OnSeq chr loc) = location t + trxlen = Loc.length loc + efflen = min inlen trxlen + in do subsubloc <- if efflen > 0 + then Just $! Loc.fromBoundsStrand (trxlen - efflen) (trxlen - 1) Plus + else Nothing + subloc <- Loc.clocOutof subsubloc loc + return $! t { location = (OnSeq chr subloc), cds = Nothing } + in test_subtranscript subtrx [ "--whole-trx", "--truncate", "--end=" ++ show (Pos.unOff outlen), "--length=" ++ show (Pos.unOff $ inlen + outlen) ] + +test_subtranscript :: (Transcript -> Maybe Transcript) -> [String] -> Property +test_subtranscript subtrx args = monadicIO $ + forAllM (listOf1 genTranscript) $ \trxs -> + let subtrxs = mapMaybe subtrx trxs + in run $ do (fullname, hfull) <- openTempFile "." "test-subtrx-full.bed" + (subname, hsub) <- openTempFile "." "test-subtrx-sub.bed" + (outname, hout) <- openTempFile "." "test-subtrx-out.bed" + hClose hout + removeFile outname + BS.hPutStr hfull . BS.unlines . map Bed.transcriptToBedStd $ trxs + hClose hfull + BS.hPutStr hsub . BS.unlines . map Bed.transcriptToBedStd $ subtrxs + hClose hsub + callProcess bedSubregion (args ++ [ "-i", fullname, "-o", outname ]) + callProcess "diff" [ subname, outname ] + removeFile fullname + removeFile subname + removeFile outname + +--test_cdsStartLength :: Property +--test_cdsStartLength = monadicIO $ +-- forAllM (listOf1 genTranscript) $ \trxs -> +-- let subtrxs = su + +-- | + +genName :: Gen SeqLabel +genName = liftM (toSeqLabel . BS.pack) $ genNameLength >>= flip replicateM genNameChar + where genNameLength = choose (1, 20) + genNameChar = elements $ ['a'..'z'] ++ ['A'..'Z'] ++ ['0'..'9'] ++ "-_" + +instance Arbitrary SeqLabel where + arbitrary = genName + +-- | Constrained position generators + +genOffset :: Gen Pos.Offset +genOffset = do isneg <- arbitrary + nnoff <- genNonNegOffset + return $ (if isneg then negate else id) nnoff + +genNonNegOffset :: Gen Pos.Offset +genNonNegOffset = liftM (subtract 1) genPositiveOffset + +genPositiveOffset :: Gen Pos.Offset +genPositiveOffset = do scale <- chooseInteger (1, 10) + liftM fromIntegral $ chooseInteger (1, 2^scale) + where chooseInteger :: (Integer, Integer) -> Gen Integer + chooseInteger = choose + +genInvertibleLoc :: Gen SpLoc.SpliceLoc +genInvertibleLoc = sized $ \sz -> do ncontigs <- choose (1, sz + 1) + fwdloc <- liftM (fromJust . SpLoc.fromContigs) + $ genContigs ncontigs + rc <- arbitrary + if rc then return $ revCompl fwdloc else return fwdloc + where genContigs = liftM (reverse . foldl' intervalsToContigs []) . genIntervals + genIntervals nints = replicateM nints $ liftM2 (,) genPositiveOffset genPositiveOffset + intervalsToContigs [] (init5, len) = [Loc.fromPosLen (Pos.Pos init5 Plus) len] + intervalsToContigs prevs@(prev:_) (nextoffset, nextlen) + = let !prevend = Loc.offset5 prev + Loc.length prev + in (Loc.fromPosLen (Pos.Pos (prevend + nextoffset) Plus) nextlen):prevs + +instance Arbitrary SpLoc.SpliceLoc where + arbitrary = genInvertibleLoc + +genTranscript :: Gen Transcript +genTranscript = do gene <- arbitrary + trx <- arbitrary + seqloc <- OnSeq <$> arbitrary <*> genInvertibleLoc + cds <- arbitrary >>= \hascds -> + if hascds + then let len = fromIntegral . Loc.length . unOnSeq $ seqloc + in do s <- Pos.Offset <$> choose (0, len - 1) + e <- Pos.Offset <$> choose (fromIntegral s + 1, len - 1) + return . Just $! Loc.fromBoundsStrand s e Plus + else return Nothing + return $! Transcript { geneId = gene, trxId = trx, location = seqloc, cds = cds } + +instance Show Transcript where + show t = "Transcript " ++ (show . BS.unpack . unSeqLabel . geneId $ t) + +data Test = forall t . Testable t => T String t + +runTest :: Test -> IO () +runTest (T name test) = do + putStr $ name ++ replicate (40 - length name) '.' ++ " " + quickCheckWith args test + where args = stdArgs -- { maxDiscard = 100000 } +