seqloc-datafiles 0.4 → 0.4.2
raw patch · 5 files changed
+341/−46 lines, 5 filesdep +directorydep +processdep ~iterateedep ~lifted-basedep ~seqloc
Dependencies added: directory, process
Dependency ranges changed: iteratee, lifted-base, seqloc, transformers-base
Files
- seqloc-datafiles.cabal +46/−27
- src/BedSubregion.hs +54/−9
- src/Bio/SeqLoc/Bed.hs +16/−4
- test/TestBed.hs +11/−6
- test/TestSubregions.hs +214/−0
seqloc-datafiles.cabal view
@@ -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
src/BedSubregion.hs view
@@ -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
src/Bio/SeqLoc/Bed.hs view
@@ -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 ])
test/TestBed.hs view
@@ -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'
+ test/TestSubregions.hs view
@@ -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 }+