twobitreader (empty) → 1.0
raw patch · 6 files changed
+1153/−0 lines, 6 filesdep +basedep +bytestringdep +directory
Dependencies added: base, bytestring, directory, mmap, optparse-applicative, primitive, twobitreader, unordered-containers, zlib
Files
- Bio/TwoBit.hs +459/−0
- Bio/TwoBit/Tool.hs +483/−0
- CHANGELOG.md +5/−0
- LICENSE +30/−0
- tool/twobit.hs +122/−0
- twobitreader.cabal +54/−0
+ Bio/TwoBit.hs view
@@ -0,0 +1,459 @@+-- | .2bit format (from the UCSC Genome Browser FAQ)+--+-- A .2bit file stores multiple DNA sequences (up to 4 Gb total) in a+-- compact randomly-accessible format. The file contains masking+-- information as well as the DNA itself.+--+-- The file begins with a 16-byte header containing the following fields:+--+-- - signature - the number 0x1A412743 in the architecture of the machine that created the file+-- - version - zero for now. Readers should abort if they see a version number higher than 0+-- - sequenceCount - the number of sequences in the file+-- - reserved - always zero for now+--+-- All fields are 32 bits unless noted. If the signature value is not as+-- given, the reader program should byte-swap the signature and check if+-- the swapped version matches. If so, all multiple-byte entities in the+-- file will have to be byte-swapped. This enables these binary files to+-- be used unchanged on different architectures.+--+-- The header is followed by a file index, which contains one entry for+-- each sequence. Each index entry contains three fields:+--+-- - nameSize - a byte containing the length of the name field+-- - name - the sequence name itself (in ASCII-compatible byte string), of variable length depending on nameSize+-- - offset - the 32-bit offset of the sequence data relative to the start of the file, not aligned to any 4-byte padding boundary+--+-- The index is followed by the sequence records, which contain nine fields:+--+-- - dnaSize - number of bases of DNA in the sequence+-- - nBlockCount - the number of blocks of Ns in the file (representing unknown sequence)+-- - nBlockStarts - an array of length nBlockCount of 32 bit integers indicating the (0-based) starting position of a block of Ns+-- - nBlockSizes - an array of length nBlockCount of 32 bit integers indicating the length of a block of Ns+-- - maskBlockCount - the number of masked (lower-case) blocks+-- - maskBlockStarts - an array of length maskBlockCount of 32 bit integers indicating the (0-based) starting position of a masked block+-- - maskBlockSizes - an array of length maskBlockCount of 32 bit integers indicating the length of a masked block+-- - reserved - always zero for now+-- - packedDna - the DNA packed to two bits per base, represented as so:+-- T - 00, C - 01, A - 10, G - 11. The first base is in the most+-- significant 2-bit byte; the last base is in the least significant+-- 2 bits. For example, the sequence TCAG is represented as 00011011.+--+-- In this format, it is neither possible nor necessary to store Ns in+-- the main sequence, and one wouldn't expect them to take up space+-- there. However, they do; hard masked sequence is typically stored as+-- many Ts. The sensible way to treat these is probably to just say+-- there are two kinds of implied annotation (repeats and large gaps for+-- a typical genome), which can be interpreted in whatever way fits.++module Bio.TwoBit (+ TwoBitFile(..),+ openTwoBit,++ TwoBitChromosome(..),+ tbf_chrnames,+ findChrom,++ TwoBitSequence'(..),+ TwoBitSequence,+ Unidrectional,+ Bidirectional,+ unpackRSRaw,+ unpackRS,+ unpackRSMasked,++ Masking(..),+ isSoftMasked,+ isHardMasked,+ noneMasked,+ softMasked,+ hardMasked,+ bothMasked+ ) where++import Control.Applicative+import Control.Exception ( Exception(..), throw )+import Control.Monad ( guard )+import Control.Monad.Primitive ( unsafeInlineIO )+import Data.Bits+import Data.Char ( toLower )+import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString.Internal as B ( fromForeignPtr )+import Data.Foldable+import qualified Data.HashMap.Strict as M+import Data.List ( mapAccumL )+import Data.Primitive.Array ( Array, arrayFromList )+import Data.Primitive.PrimArray ( indexPrimArray )+import Data.Word ( byteSwap32, Word8, Word32 )+import Foreign.ForeignPtr ( ForeignPtr, withForeignPtr )+import Foreign.Ptr ( castPtr, plusPtr, Ptr )+import Foreign.Storable ( Storable(..) )+import GHC.Base ( build )+import System.IO.MMap ( mmapFileForeignPtr, Mode(..) )++data TwoBitFile = TBF { tbf_raw :: {-# UNPACK #-} !(ForeignPtr Word8)+ , tbf_size :: {-# UNPACK #-} !Int+ , tbf_path :: {-# UNPACK #-} !B.ByteString+ , tbf_chroms :: {-# UNPACK #-} !(Array TwoBitChromosome)+ , tbf_chrmap :: !(M.HashMap B.ByteString TwoBitChromosome) }++tbf_chrnames :: TwoBitFile -> [B.ByteString]+tbf_chrnames = toList . fmap tbc_name . tbf_chroms++-- | Finds a named scaffold in the reference. If it doesn't find the+-- exact name, it will try to compensate for the crazy naming+-- differences between NCBI and UCSC. This doesn't work in general, but+-- is good enough in the common case. In particular, "1" maps to "chr1"+-- and back, "GL000192.1" to "chr1_gl000192_random" and back, and "chrM"+-- to "MT" and back.+findChrom :: B.ByteString -> TwoBitFile -> Maybe TwoBitChromosome+findChrom c TBF{ tbf_chrmap = cs } =+ M.lookup c cs+ <|> M.lookup ("chr" <> c) cs+ <|> ( guard ("chr" `B.isPrefixOf` c) >> M.lookup (B.drop 3 c) cs )+ <|> ( guard ("chrM" == c) >> M.lookup "MT" cs )+ <|> ( guard ("MT" == c) >> M.lookup "chrM" cs )+ <|> ( case filter (\d -> match c (tbc_name d) || match (tbc_name d) c) $ M.elems cs of+ [x] -> Just x ; _ -> Nothing )+ where+ match x y = B.isInfixOf (B.map toLower (B.takeWhile (/= '.') x)) y++data TwoBitChromosome = TBC { tbc_raw :: {-# UNPACK #-} !(ForeignPtr Word8)+ , tbc_name :: {-# UNPACK #-} !B.ByteString+ , tbc_index :: {-# UNPACK #-} !Int+ , tbc_dna_offset :: {-# UNPACK #-} !Word32+ , tbc_dna_size :: {-# UNPACK #-} !Word32+ -- | Lazily generated sequence in forward direction; the argument is the offset of the first base.+ , tbc_fwd_seq :: Int -> TwoBitSequence' Unidrectional+ -- | Lazily generated sequence in reverse direction; the argument is the offset of the first base to the+ -- right of the beginning. (The first base generated is the complement of the base found at (offset-1).+ , tbc_rev_seq :: Int -> TwoBitSequence' Bidirectional }+++data TwoBitError = WrongSignature FilePath+ | UnsortedBlocks FilePath+ | OutOfBounds FilePath Word32 Int+ | OverlongSequence FilePath Word32 Word32 Int+ deriving Show++instance Exception TwoBitError where+ displayException (WrongSignature fp) = "The file " ++ show fp ++ "does not have a .2bit signature."+ displayException (UnsortedBlocks fp) = "The N and mask blocks in file " ++ show fp ++ " are not sorted."+ displayException (OutOfBounds fp o s) = "Attempted to access offset " ++ show o ++ " in file " ++ show fp ++ " of size " ++ show s ++ "."+ displayException (OverlongSequence fp o l s) = "A sequence of length " ++ show l ++ " starting at " ++ show o ++ " in file "+ ++ show fp ++ " hangs over its end at " ++ show s ++ "."+++-- | Brings a 2bit file into memory. The file is mmap'ed, so it will+-- not work on streams that are not actual files. It's also unsafe if+-- the file is concurrently modified in any way.+openTwoBit :: FilePath -> IO TwoBitFile+openTwoBit fp = do (p,o,l) <- mmapFileForeignPtr fp ReadOnly Nothing+ if o == 0 then pure $ parseTwoBit fp p l+ else fail $ "unexpected: mmapFileForeignPtr returned an offset"+++-- | Parses a 2bit file. The @FilePath@ argument is only used in error+-- messages, what is really parsed is the memory block, typically from+-- mmapping the file.+--+-- The workhorse in here is the construction of the 'tbc_fwd_seq' and+-- 'tbc_rev_seq' functions. When called, they first run a binary search+-- on the mask lists, then produce a list of blocks with uniform+-- masking. Both parts of the algorithm are fast and directly use the+-- on-disk data structures.+--+-- In theory, there could be 2bit files in big endian format out there.+-- We nominally support them, but since I've never seen one in the wild,+-- this may well fail in a spectacular way.++parseTwoBit :: FilePath -> ForeignPtr Word8 -> Int -> TwoBitFile+parseTwoBit fp0 raw size+ | getW32_ peekUnalnWord32 0 == 0x1A412743 && getW32_ peekUnalnWord32 4 == 0 = kont $ parseEachSeq (getW32_ peekUnalnWord32)+ | getW32_ peekUnalnWord32Swap 0 == 0x1A412743 && getW32_ peekUnalnWord32Swap 4 == 0 = kont $ parseEachSeq (getW32_ peekUnalnWord32Swap)+ | otherwise = throw $ WrongSignature fp0+ where+ kont sqs = TBF raw size (B.pack fp0) (arrayFromList sqs) (M.fromList $ map (liftA2 (,) tbc_name id) sqs)++ getW32_ f o | o + 4 >= fromIntegral size = throw $ OutOfBounds fp0 o size+ | otherwise = unsafeInlineIO $ withForeignPtr raw $ \p -> f (plusPtr p (fromIntegral o))++ parseEachSeq :: (Word32 -> Word32) -> [TwoBitChromosome]+ parseEachSeq getW32 = snd $ mapAccumL (parseOneSeq getW32) 16 [0 .. fromIntegral (getW32 8) -1]++ parseOneSeq getW32 off nseq =+ if packedDnaOff + shiftR (dnasize+3) 2 > fromIntegral size+ then throw $ OverlongSequence fp0 packedDnaOff dnasize size+ else (off + 5 + nmsize, TBC raw name nseq packedDnaOff dnasize unfoldSeqFwd unfoldSeqRev)+ where+ !nmsize = unsafeInlineIO $ withForeignPtr raw $ \p -> fromIntegral <$> peekElemOff p off+ !name = B.fromForeignPtr raw (off+1) nmsize+ !offset = getW32 . fromIntegral $ off + 1 + nmsize++ !dnasize = getW32 $ offset+ !nBlockCount = getW32 $ offset + 4+ !mBlockCount = getW32 $ offset + 8 + 8*nBlockCount+ !packedDnaOff = offset + 16 + 8 * (nBlockCount+mBlockCount)++ -- Valid blocks are numbered 1..max; there are virtual guard blocks at indices 0 and (max+1), which make the later+ -- algorithms much cleaner+ n_block, m_block :: Word32 -> Block+ n_block i | i == 0 = B 0 0 i+ | i > nBlockCount = B maxBound maxBound i+ | otherwise = B a (a+b) i+ where+ !a = getW32 $ offset+4 + 4*i+ !b = getW32 $ offset+4 + 4*(i+nBlockCount)++ m_block i | i == 0 = B 0 0 i+ | i > mBlockCount = B maxBound maxBound i+ | otherwise = B a (a+b) i+ where+ !a = getW32 $ offset+8 + 8*nBlockCount + 4*i+ !b = getW32 $ offset+8 + 8*nBlockCount + 4*(i+mBlockCount)+++ unfoldSeqFwd :: Int -> TwoBitSequence' Unidrectional+ unfoldSeqFwd chroff = unfoldSeqFwd' (search n_block nBlockCount) (search m_block mBlockCount)+ (fromIntegral chroff) (packedDnaOff * 4 + fromIntegral chroff)+ where+ trim b = b { start_offset = max (fromIntegral chroff) (start_offset b) }++ -- finds the smallest index such that the block end(!) is larger than 'chroff'+ search f num = trim . f $ go 0 (num+1)+ where+ go a b | a == b = a+ | end_offset (f m) > fromIntegral chroff = go a m+ | otherwise = go (m+1) b+ where+ m = div (a + b) 2++ unfoldSeqFwd' :: Block -> Block -> Word32 -> Word32 -> TwoBitSequence' Unidrectional+ unfoldSeqFwd' nb@(B nstart nend _) mb@(B mstart mend _) !chroff !fileoff+ | chroff >= dnasize = RefEnd+ | chroff > mstart || chroff > nstart = throw (UnsortedBlocks fp0)+ | chroff < nstart && chroff < mstart = advance noneMasked $ min dnasize $ min nstart mstart+ | chroff < mstart = advance hardMasked $ min dnasize $ min nend mstart+ | chroff < nstart = advance softMasked $ min dnasize $ min mend nstart+ | otherwise = advance bothMasked $ min dnasize $ min nend mend+ where+ advance m x = SomeSeq m raw (fromIntegral fileoff) (fromIntegral $ x - chroff) $+ unfoldSeqFwd' (advanceB x n_block nb) (advanceB x m_block mb) x (fileoff+x-chroff)++ advanceB :: Word32 -> (Word32 -> Block) -> Block -> Block+ advanceB x f (B start end i)+ | x <= start = B start end i+ | x < end = B x end i+ | otherwise = f (i+1)+++ unfoldSeqRev :: Int -> TwoBitSequence' Bidirectional+ unfoldSeqRev chroff = unfoldSeqRev' (search n_block nBlockCount) (search m_block mBlockCount)+ (fromIntegral chroff) (packedDnaOff * 4 + fromIntegral chroff)+ where+ trim b = b { end_offset = min (fromIntegral chroff) (end_offset b) }++ -- finds the largest index such that the block start is smaller than chroff+ search f num = trim . f $ go 0 (num+1)+ where+ go a b | a == b = a+ | start_offset (f m) < fromIntegral chroff = go m b+ | otherwise = go a (m-1)+ where+ m = div (a + b + 1) 2++ unfoldSeqRev' :: Block -> Block -> Word32 -> Word32 -> TwoBitSequence' Bidirectional+ unfoldSeqRev' nb@(B nstart nend _) mb@(B mstart mend _) !chroff !fileoff+ | chroff <= 0 = RefEnd+ | chroff < mend || chroff < nend = throw (UnsortedBlocks fp0)+ | chroff > nend && chroff > mend = advance noneMasked $ max nend mend+ | chroff > mend = advance hardMasked $ max nstart mend+ | chroff > nend = advance softMasked $ max mstart nend+ | otherwise = advance bothMasked $ max nstart mstart+ where+ advance m x = SomeSeq m raw (fromIntegral fileoff) (fromIntegral x - fromIntegral chroff) $+ unfoldSeqRev' (advanceB x n_block nb) (advanceB x m_block mb) x (fileoff+x-chroff)++ advanceB :: Word32 -> (Word32 -> Block) -> Block -> Block+ advanceB x f (B start end i)+ | x >= end = B start end i+ | x > start = B start x i+ | otherwise = f (i-1)+++data Block = B { start_offset :: !Word32+ , end_offset :: !Word32+ , block_number :: !Word32 }+ deriving (Show, Eq, Ord)+++-- | 2bit supports two kinds of masking, typically rendered as lowercase+-- letters ('MaskSoft') and Ns ('MaskHard'). They can overlap+-- ('MaskBoth'), and even the hard masking has underlying sequence+-- (which is normally ignored).+newtype Masking = Masking Word8 deriving (Eq, Ord)++instance Show Masking where+ show (Masking 0) = "None"+ show (Masking 1) = "Soft"+ show (Masking 2) = "Hard"+ show (Masking _) = "Both"++instance Read Masking where+ readsPrec _ s = [ (Masking m,s') | (w,s') <- lex s+ , m <- case w of "None" -> [0]+ "Soft" -> [1]+ "Hard" -> [2]+ "Both" -> [3]+ _ -> [ ] ]++instance Semigroup Masking where+ Masking a <> Masking b = Masking (a .|. b)++instance Monoid Masking where+ mempty = Masking 0+ mappend = (<>)++instance Enum Masking where+ toEnum = Masking . toEnum+ fromEnum (Masking m) = fromEnum m++instance Bounded Masking where+ minBound = Masking 0+ maxBound = Masking 3+++isSoftMasked, isHardMasked :: Masking -> Bool+isSoftMasked (Masking m) = testBit m 0+isHardMasked (Masking m) = testBit m 1++noneMasked, softMasked, hardMasked, bothMasked :: Masking+noneMasked = Masking 0+softMasked = Masking 1+hardMasked = Masking 2+bothMasked = Masking 3+++-- | This is a (piece of a) reference sequence. It consists of+-- stretches with uniform masking.+--+-- The offset is stored as a 'Word'. This is done because on a 32 bit+-- platform, every bit counts. This limits the genome to approximately+-- four gigabases, which would be a file of about one gigabyte. That's+-- just about enough to work with the human genome. On a 64 bit+-- platform, the file format itself imposes a limit of four gigabytes,+-- or about 16 gigabases in total.+--+-- If length is zero, the piece is empty and the mask, pointer, and+-- offset fields may not be valid. If length is positive, ptr+offset+-- points at the first base of the piece. If length is negative,+-- ptr+offset points just past the end of the piece, ptr+offset+length+-- points to the first base of the piece, and the sequence in meant to+-- be reverse complemented.+--+-- In a 'TwoBitSequence', length must not be negative. In a+-- @TwoBitSequence' Bidirectional@, length can be positive or negative.++data TwoBitSequence' dir = SomeSeq {-# UNPACK #-} !Masking -- ^ how is it masked?+ {-# UNPACK #-} !(ForeignPtr Word8) -- ^ primitive bases in 2bit encoding: [0..3] = TCAG+ {-# UNPACK #-} !Word -- ^ offset in bases(!)+ {-# UNPACK #-} !Int -- ^ length in bases+ (TwoBitSequence' dir)+ | RefEnd++data Unidrectional+data Bidirectional++type TwoBitSequence = TwoBitSequence' Unidrectional++instance Show (TwoBitSequence' dir) where+ showsPrec _ (SomeSeq m _ _ l r) = (++) "SomeSeq " . shows m . (:) ' ' . shows l . (++) " $ " . shows r+ showsPrec _ RefEnd = (++) "RefEnd"++-- | Unpacks a reference sequence into a (very long) list of bytes.+-- Each byte contains the nucleotide in bits 0 and 1 with valjues 0..3+-- corresponding to "TCAG", and the soft and hard mask bits in bits 2+-- and 3, respectively.++unpackRSRaw :: TwoBitSequence' dir -> [Word8]+unpackRSRaw rs = build (\c n -> unpackRSFB c n rs)+{-# INLINE unpackRSRaw #-}++unpackRSFB :: (Word8 -> b -> b) -> b -> TwoBitSequence' dir -> b+unpackRSFB cons nil = go0+ where+ go0 RefEnd = nil+ go0 (SomeSeq (Masking msk) raw off0 len0 rs) | len0 >= 0 = go off0 len0+ where+ go !off !len = if len == 0 then go0 rs else code `cons` go (off+1) (len-1)+ where+ !byteoff = fromIntegral $ off `shiftR` 2+ !bitoff = fromIntegral $ off .&. 3+ !byte = unsafeInlineIO $ withForeignPtr raw (`peekByteOff` byteoff)+ !code = shiftR byte (6 - 2 * bitoff ) .&. 3 .|. shiftL msk 2++ go0 (SomeSeq (Masking msk) raw off0 len0 rs) = go off0 (-len0)+ where+ go !off !len = if len == 0 then go0 rs else xor 2 code `cons` go (off-1) (len-1)+ where+ !byteoff = fromIntegral $ (off-1) `shiftR` 2+ !bitoff = fromIntegral $ (off-1) .&. 3+ !byte = unsafeInlineIO $ withForeignPtr raw (`peekByteOff` byteoff)+ !code = shiftR byte (6 - 2 * bitoff ) .&. 3 .|. shiftL msk 2+{-# INLINE [0] unpackRSFB #-}++-- | Unpacks a reference sequence into a (very long) list of ASCII+-- characters. Hard masked nucleotides become the letter 'N', others+-- become "TCAG".+unpackRS :: TwoBitSequence' dir -> [Word8]+unpackRS = map (indexPrimArray chars . fromIntegral) . unpackRSRaw+ where+ !chars = [84,67,65,71,84,67,65,71,78,78,78,78,78,78,78,78] -- "TCAGTCAGNNNNNNNN"+{-# INLINE unpackRS #-}++-- | Unpacks a reference sequence into a list of ASCII characters,+-- interpreting masking in the customary way. Specifically, hard+-- masking produces Ns, soft masking produces lower case letters, and+-- dual masking produces lower case Ns.+unpackRSMasked :: TwoBitSequence' dir -> [Word8]+unpackRSMasked = map (indexPrimArray chars . fromIntegral) . unpackRSRaw+ where+ !chars = [84,67,65,71,116,99,97,103,78,78,78,78,110,110,110,110] -- "TCAGtcagNNNNnnnn"+{-# INLINE unpackRSMasked #-}+++-- | Reads a 32 bit word from an address, which doesn't need to be+-- aligned. The byte order used is unspecified.+peekUnalnWord32 :: Ptr a -> IO Word32++-- | Equivalent to peekUnalnWord32 followed by a byte swap.+peekUnalnWord32Swap :: Ptr a -> IO Word32+++-- List of known architectures that efficiently support unaligned accesses.+#if defined(i386_HOST_ARCH) || defined(x86_64_HOST_ARCH) \+ || defined(powerpc64le_HOST_ARCH) || ((defined(arm_HOST_ARCH) \+ || defined(aarch64_HOST_ARCH)) && defined(__ARM_FEATURE_UNALIGNED)) \+ || defined(powerpc_HOST_ARCH) || defined(powerpc64_HOST_ARCH)++peekUnalnWord32 = peek . castPtr+peekUnalnWord32Swap = fmap byteSwap32 . peek . castPtr++#else++peekUnalnWord32 p = do+ x <- fromIntegral <$> peekWord8 (plusPtr p 0)+ y <- fromIntegral <$> peekWord8 (plusPtr p 1)+ z <- fromIntegral <$> peekWord8 (plusPtr p 2)+ w <- fromIntegral <$> peekWord8 (plusPtr p 3)+ return $! x .|. unsafeShiftL y 8 .|. unsafeShiftL z 16 .|. unsafeShiftL w 24++peekUnalnWord32Swap p = do+ x <- fromIntegral <$> peekWord8 (plusPtr p 0)+ y <- fromIntegral <$> peekWord8 (plusPtr p 1)+ z <- fromIntegral <$> peekWord8 (plusPtr p 2)+ w <- fromIntegral <$> peekWord8 (plusPtr p 3)+ return $! w .|. unsafeShiftL z 8 .|. unsafeShiftL y 16 .|. unsafeShiftL x 24++#endif+
+ Bio/TwoBit/Tool.hs view
@@ -0,0 +1,483 @@+{-# OPTIONS_GHC -Wno-partial-fields #-}+module Bio.TwoBit.Tool+ ( EncodeProgress(..)+ , buildFasta+ , faToTwoBit+ , formatCdna+ , parseAnno+ , twoBitToFa+ , vcfToTwoBit+ )+where++import Bio.TwoBit+import Control.Applicative+import Control.Exception+import Control.Monad+import Data.Bits+import Data.Bool+import qualified Data.ByteString as B+import qualified Data.ByteString.Builder as B+import qualified Data.ByteString.Char8 as C+import qualified Data.ByteString.Lazy.Char8 as L+import Data.ByteString.Short ( ShortByteString, toShort )+import qualified Data.ByteString.Short as H+import Data.Char ( isSpace, isUpper )+import Data.Foldable+import qualified Data.HashMap.Strict as M+import Data.Int ( Int64 )+import Data.Word ( Word8, Word32 )+import System.IO ( stdout )++type Bytes = B.ByteString+type LazyBytes = L.ByteString++-- | A cDNA or mRNA or transcript (these are all synonymous), with some+-- metainformation collected from the annotation. Whatever the input+-- was called, we call it 'cdna' in the transciptome.+data Cdna = Cdna+ { c_id :: !Bytes -- identifier (typically an ENST number)+ , c_pos :: !Range -- genomic position+ , c_gene_id :: !Bytes -- gene identifier (typically an ENSG number)+ , c_gene_symbol :: !Bytes -- colloquial name, aka locus+ , c_gene_biotype :: !Bytes -- whatever, just pass it on+ , c_biotype :: !Bytes -- whatever, just pass it on+ , c_description :: !Bytes -- unclear; always empty for now+ , c_exons :: [Range] -- list of exon coordinates (sorted backwards)+ }+ deriving Show++data Range = Range+ { r_chrom :: !C.ByteString+ , r_start :: !Int+ , r_len :: !Int }+ deriving Show++reverseRange :: Range -> Range+reverseRange (Range sq pos len) = Range sq (-pos-len) len++null_cdna :: Cdna+null_cdna = Cdna "" (Range "" 0 0) "" "" "" "" "" []+++formatCdna :: TwoBitFile -> Cdna -> B.Builder+formatCdna tbf Cdna{..} = descr <> buildFasta 60 getExons+ where+ (_,tbf_fn) = C.breakEnd (=='/') $ tbf_path tbf+ (tbf_base,_) = C.breakEnd (=='.') tbf_fn+ genome_id = if C.null tbf_base then tbf_fn else C.init tbf_base++ descr = B.char7 '>' <> B.byteString c_id <> " cdna chromosome:" <>+ B.byteString genome_id <> B.char7 ':' <> formatRange c_pos <>+ " gene:" <> B.byteString c_gene_id <>+ maybeBS " gene_biotype:" c_gene_biotype <>+ maybeBS " transcript_biotype:" c_biotype <>+ maybeBS " gene_symbol:" c_gene_symbol <>+ maybeBS " description:" c_description <> B.char7 '\n'++ formatRange r | r_start r < 0 = formatRange1 (reverseRange r) <> ":-1"+ | otherwise = formatRange1 r <> ":1"++ formatRange1 r = B.byteString (r_chrom r) <> B.char7 ':' <>+ B.intDec (r_start r) <> B.char7 ':' <>+ B.intDec (r_start r + r_len r - 1)++ maybeBS p s = if B.null s then mempty else p <> B.byteString s++ getExons | r_start c_pos < 0 = concatMap getExon c_exons+ | otherwise = concatMap getExon (reverse c_exons)++ getExon :: Range -> [Word8]+ getExon (Range ch start len) =+ case findChrom ch tbf of+ Just tbs | start >= 0 -> take len $ unpackRS $ tbc_fwd_seq tbs start+ | otherwise -> take len $ unpackRS $ tbc_rev_seq tbs (-start-len)+ Nothing -> error $ "unknown reference " ++ show ch+++++buildFasta :: Int -> [Word8] -> B.Builder+buildFasta n = go+ where+ go [ ] = mempty+ go s = let (u,v) = splitAt n s+ in foldMap B.word8 u <> B.char7 '\n' <> go v+{-# INLINE buildFasta #-}++twoBitToFa :: Int -> TwoBitSequence' dir -> IO ()+twoBitToFa ln = B.hPutBuilder stdout . buildFasta 60 . take ln . unpackRSMasked+++data EncodeProgress+ = EncodeProgress+ { ep_seqname :: !ShortByteString+ , ep_position :: !Word32+ , ep_hardmasked :: !Word32+ , ep_softmasked :: !Word32+ , ep_enclength :: !Int64+ , ep_tail :: EncodeProgress }+ | Encoded B.Builder+++-- Strategy: We can only write the packedDNA after we wrote the nBlocks+-- and mBlocks. So packedDNA needs to be buffered. We have to do three+-- simultaneous strict folds of the input, all of which result in reasonably+-- compact structures (name table, mask table, encoded dna), which get+-- concatenated at the end.+--+-- We also have to buffer everything, since the header with the sequence+-- names must be written first. Oh joy.+--+-- We return a list of progress notifications terminated by the+-- 'B.Builder' for the whole 2bit file. The progress messages can be+-- printed or ignored; in either case, they should ensure enough+-- strictness to not waste more memory than necessary.++faToTwoBit :: L.ByteString -> EncodeProgress+faToTwoBit = get_each []+ where+ get_each acc inp = case L.uncons $ L.dropWhile (/= '>') inp of+ Nothing -> Encoded $ seqs_to_twobit $ reverse acc+ Just (_,s2) ->+ let (nm, s') = L.break (<= ' ') s2+ in get_one acc (toShort (L.toStrict nm)) 0 (GapList maxBound L2i_Nil)+ (GapList maxBound L2i_Nil) (BaseAccu 0 0 emptyAccu)+ (L.dropWhile (/= '\n') s')++ get_one acc !nm !pos !ns !ms !bs inp = case L.uncons inp of+ Nothing -> fin L.empty+ Just (c,s')+ | c <= ' ' -> get_one acc nm pos ns ms bs s'+ | c == '>' -> fin (L.cons c s')+ | otherwise -> get_one acc nm (succ pos)+ (collect_Ns ns pos c)+ (collect_ms ms pos c)+ (collect_bases bs c) s'+ where+ fin k = let !r = encode_seq pos ns ms bs+ in EncodeProgress nm pos (sum_L2i pos ns) (sum_L2i pos ms) (L.length r) $+ get_each ((nm,r):acc) k++-- | Extracts the reference from a VCF. This assumes the presence of at+-- least one record per site. The VCF must be sorted by position. When+-- writing out, we try to match the order of the contigs as listed in+-- the header. Unlisted contigs follow at the end with their order+-- preserved; contigs without data are not written at all.+vcfToTwoBit :: [B.ByteString] -> EncodeProgress+vcfToTwoBit s0 = let (lns, s1) = read_header [] s0+ in get_each lns [] $ filter (\s -> not (B.null s) && C.head s /= '#') s1+ where+ -- Collects the "contig" stanzas, parses their lengths. Returns the+ -- length map and the remaining stream.+ read_header acc [ ] = (reverse acc, [])+ read_header acc (l:ls) | "##contig=" `C.isPrefixOf` l+ , (Just !nm, Just !ln) <- parse_cline l = read_header ((nm,ln):acc) ls+ | "#" `C.isPrefixOf` l = read_header acc ls+ | otherwise = (reverse acc, l:ls)++ parse_cline = p1 . C.filter (not . isSpace) . C.takeWhile (/='>') . C.drop 1 . C.dropWhile (/='<')+ where+ p1 s | "ID=" `C.isPrefixOf` s = let (nm,t) = C.break (==',') $ C.drop 3 s+ (_,ln) = p1 $ C.drop 1 t+ in (Just (toShort nm),ln)++ | "length=" `C.isPrefixOf` s = case C.readInt $ C.drop 7 s of+ Just (ln,u) -> let (nm,_) = p1 $ C.drop 1 $ C.dropWhile (/=',') u in (nm,Just (fromIntegral ln))+ Nothing -> p1 $ C.drop 1 $ C.dropWhile (/=',') s++ | C.null s = (Nothing,Nothing)+ | otherwise = p1 $ C.drop 1 $ C.dropWhile (/=',') s++ get_each :: [(ShortByteString,Word32)]+ -> [(ShortByteString, LazyBytes)]+ -> [B.ByteString]+ -> EncodeProgress+ get_each lns acc [ ] = Encoded $ seqs_to_twobit $ reorder (map fst lns) $ reverse acc+ get_each lns acc (l:s2) = EncodeProgress nm' ln' (sum_L2i ln' ns') 0 (L.length r) $+ get_each lns ((nm',r):acc) s3+ where+ nm = B.takeWhile (/= 9) l+ (pos,ns,bs,s3) = get_one nm 0 (GapList maxBound L2i_Nil) (BaseAccu 0 0 emptyAccu) (l:s2)+ !nm' = toShort nm+ (ns',bs',ln') = case find ((==) nm' . fst) lns of+ Just (_,ln) | ln > pos -> (extend_gap ns ln, pad_bases bs (fromIntegral $ ln-pos), ln)+ _ -> (ns,bs,pos)+ !r = encode_seq ln' ns' (GapList maxBound L2i_Nil) bs'+++ -- important: 1-based coordinates!+ get_one !_nm !pos !ns !bs [ ] = (pos,ns,bs,[])+ get_one !nm !pos !ns !bs (l:s')+ | B.takeWhile (/=9) l /= nm = (pos,ns,bs,l:s')++ | Just (pos',l3) <- C.readInt . C.drop 1 $ B.dropWhile (/=9) l+ , ref <- B.takeWhile (/=9) . B.drop 1 . B.dropWhile (/=9) $ B.drop 1 l3+ , fromIntegral pos' >= pos + 1+ , not (C.null ref) =+ if fromIntegral pos' == pos + 1+ -- record in sequence+ then get_one nm (succ pos) (collect_Ns ns pos $ C.head ref)+ (collect_bases bs $ C.head ref) s'+ -- gap: handle the gap, reprocess the record+ else let gap_len = pos' - fromIntegral pos - 1+ in get_one nm (fromIntegral pos' - 1) (extend_gap ns pos)+ (pad_bases bs gap_len) (l:s')++ -- anything else can be ignored (parse errors or additional records)+ | otherwise = get_one nm pos ns bs s'+++ pad_bases bs n = foldl' collect_bases bs $ replicate n 'T'++ -- Reorder a key-value list so it matches the order of a list of+ -- keys. Missing keys are ignored, leftover pairs retain their+ -- original order.+ reorder :: Eq a => [a] -> [(a,b)] -> [(a,b)]+ reorder [ ] vs = vs+ reorder (k:ks) vs = go [] vs+ where+ go xs ((k1,v1):ys) | k == k1 = (k1,v1) : reorder ks (reverse xs ++ ys)+ | otherwise = go ((k1,v1):xs) ys+ go xs [ ] = reorder ks (reverse xs)+++-- List of pairs of 'Word32's. Specialized and unpacked to conserve space. Probably overkill...+data L2i = L2i {-# UNPACK #-} !Word32 {-# UNPACK #-} !Word32 L2i | L2i_Nil++data GapList = GapList !Word32 !L2i++sum_L2i :: Word32 -> GapList -> Word32+sum_L2i p (GapList q xs) = go (if q == maxBound then 0 else p-q) xs+ where+ go !a (L2i x y z) = go (a+y-x) z+ go !a L2i_Nil = a++encodeL2i :: L2i -> B.Builder+encodeL2i = go 0 mempty mempty+ where+ go !n ss ls L2i_Nil = B.word32LE n <> ss <> ls+ go !n ss ls (L2i s e rs) = go (succ n) (B.word32LE s <> ss) (B.word32LE (e-s) <> ls) rs++seqs_to_twobit :: [(ShortByteString, LazyBytes)] -> B.Builder+seqs_to_twobit seqs = B.word32LE 0x1A412743 <> B.word32LE 0 <>+ B.word32LE (fromIntegral $ length seqs) <> B.word32LE 0 <>+ mconcat (zipWith (\nm off -> B.word8 (fromIntegral (H.length nm)) <>+ B.shortByteString nm <>+ B.word32LE (fromIntegral off))+ (map fst seqs) offsets) <>+ foldMap (B.lazyByteString . snd) seqs+ where+ offset0 = 16 + 5 * length seqs + sum (map (H.length . fst) seqs)+ offsets = scanl (\a b -> a + fromIntegral (L.length b)) offset0 $ map snd seqs+++-- | A way to accumulate bytes. If the accumulated bytes will hang+-- around in memory, this has much lower overhead than 'Builder'. If it+-- has short lifetime, 'Builder' is much more convenient.+newtype Accu = Accu [Bytes]++emptyAccu :: Accu+emptyAccu = Accu []++-- | Appends bytes to a collection of 'Bytes' in such a way that the+-- 'Bytes' keep doubling in size. This ensures O(n) time and space+-- complexity and fairly low overhead.+grow :: Word8 -> Accu -> Accu+grow w = go 1 [B.singleton w]+ where+ go l acc (Accu (s:ss))+ | B.length s <= l = go (l+B.length s) (s:acc) (Accu ss)+ | otherwise = let !s' = B.concat acc in Accu (s' : s : ss)+ go _ acc (Accu [ ]) = let !s' = B.concat acc in Accu [s']++buildAccu :: Accu -> B.Builder+buildAccu (Accu ss) = foldMap B.byteString $ reverse ss++encode_seq :: Word32 -- ^ length+ -> GapList -- ^ list of N stretches+ -> GapList -- ^ list of mask stretches+ -> BaseAccu -- ^ accumulated bases+ -> LazyBytes++encode_seq pos ns ms bs = L.length r `seq` r+ where+ ss' = case bs of (BaseAccu 0 _ ss) -> ss+ (BaseAccu n w ss) -> grow (w `shiftL` (8-2*n)) ss+ r = B.toLazyByteString $+ B.word32LE pos <>+ encodeL2i (case ns of GapList p rs | p == maxBound -> rs ; GapList p rs -> L2i p pos rs) <>+ encodeL2i (case ms of GapList p rs | p == maxBound -> rs ; GapList p rs -> L2i p pos rs) <>+ B.word32LE 0 <>+ buildAccu ss'++-- | Collects stretches of Ns by looking at one character at a time. In+-- reality, anything that isn't one of \"ACGT\" is treated as an N.+collect_Ns :: GapList -> Word32 -> Char -> GapList+collect_Ns (GapList spos rs) pos c+ | spos == maxBound && c `C.elem` "ACGTacgt" = GapList maxBound rs+ | spos == maxBound = GapList pos rs+ | c `C.elem` "ACGTacgt" = GapList maxBound (L2i spos pos rs)+ | otherwise = GapList spos rs++-- | Collects stretches of masked dna by looking at one letter at a+-- time. Anything lowercase is considered masked.+collect_ms :: GapList -> Word32 -> Char -> GapList+collect_ms (GapList spos rs) pos c+ | spos == maxBound && isUpper c = GapList maxBound rs+ | spos == maxBound = GapList pos rs+ | isUpper c = GapList maxBound (L2i spos pos rs)+ | otherwise = GapList spos rs++extend_gap :: GapList -> Word32 -> GapList+extend_gap (GapList spos rs) pos+ | spos == maxBound = GapList pos rs+ | otherwise = GapList spos rs+++data BaseAccu = BaseAccu !Int !Word8 !Accu++-- | Collects bases in 2bit format. It accumulates 4 bases in one word,+-- then collects bytes in an 'Accu'. From the 2bit spec:+--+-- packedDna - the DNA packed to two bits per base, represented as+-- so: T - 00, C - 01, A - 10, G - 11. The first base is+-- in the most significant 2-bit byte; the last base is+-- in the least significant 2 bits. For example, the+-- sequence TCAG is represented as 00011011.+collect_bases :: BaseAccu -> Char -> BaseAccu+collect_bases (BaseAccu n w ss) c+ = let code = case c of 'C'->1;'c'->1;'A'->2;'a'->2;'G'->3;'g'->3;_->0+ w' = shiftL w 2 .|. code+ in if n == 3 then BaseAccu 0 0 (grow w' ss) else BaseAccu (succ n) w' ss++data Gene = Gene { g_id :: Bytes, g_symbol :: Bytes, g_biotype :: Bytes }++null_gene :: Gene+null_gene = Gene "" "" ""++data GffError = GffError String Int GffErrorDetail deriving Show+data GffErrorDetail = GffParseError | GffIdMismatch | GffUnknownRef Bytes deriving Show++instance Exception GffError where+ displayException (GffError fp ln dt) = displayDetail dt ++ " in line " ++ show ln ++ " of gff file " ++ fp+ where+ displayDetail GffParseError = "parse error"+ displayDetail GffIdMismatch = "identifier does not match"+ displayDetail (GffUnknownRef ch) = "unknown reference " ++ show ch++-- | Parses annotations in GFF format. We want to turn an annotation+-- and a 2bit file into a FastA of the transcriptome (one sequence per+-- annotated transcript), that looks like the stuff Lior Pachter feeds+-- into Kallisto. Annotations come in two dialects of GFF, either GFF3+-- or GTF. We autodetect and understand both.++parseAnno :: String -> L.ByteString -> [Either GffError Cdna]+parseAnno fp = filter (either (const True) (not . null . c_exons)) .+ go null_gene null_cdna .+ map (fmap (C.split '\t')) .+ filter (\(_,s) -> not (B.null s) && C.head s /= '#') .+ zip (enumFrom 1) .+ map L.toStrict .+ L.lines+ where+ go gene xscript ((ln, ch:_:tp:fro_:tho_:_:strand:_:stuff_:_) : strm)+ | Just (fro,"") <- C.readInt fro_+ , Just (tho,"") <- C.readInt tho_+ , Just stuff <- parseStuff stuff_ =+ let rng = bool id reverseRange (strand == "-") $ Range ch (fro-1) (tho-fro+1)+ in go2 ln gene xscript strm rng (B.map (.|. 32) tp) stuff++ go gene xscript ((ln, _) : strm) = Left (GffError fp ln GffParseError) : go gene xscript strm+ go _gene xscript [ ] = Right xscript : []++ go2 ln gene xscript strm rng tp (Left stuff)+ | tp == "exon" = if M.lookup "Parent" stuff == Just (c_id xscript)+ then go gene xscript { c_exons = rng : c_exons xscript } strm+ else Left (GffError fp ln GffIdMismatch) : go gene xscript strm++ | tp == "transcript" || tp == "cdna" || tp == "mrna" =+ Right xscript :+ case (M.lookup "ID" stuff, M.lookup "Parent" stuff) of+ (Just tid, Just gid)+ | gid == g_id gene -> let xscript' = Cdna { c_id = tid+ , c_pos = rng+ , c_gene_id = gid+ , c_gene_symbol = g_symbol gene+ , c_gene_biotype = g_biotype gene+ , c_biotype = M.lookupDefault "" "biotype" stuff+ , c_description = "" -- XXX+ , c_exons = [] }+ in go gene xscript' strm++ | otherwise -> Left (GffError fp ln GffIdMismatch) : go gene null_cdna strm++ _ -> Left (GffError fp ln GffParseError) : go gene null_cdna strm++ | tp == "gene" =+ Right xscript :+ case M.lookup "ID" stuff of+ Just gid -> let gene' = Gene { g_id = gid+ , g_symbol = "" -- XXX+ , g_biotype = M.lookupDefault "" "biotype" stuff }+ in go gene' null_cdna strm++ Nothing -> Left (GffError fp ln GffParseError) : go null_gene null_cdna strm++ | otherwise = go gene xscript strm++ go2 ln gene xscript strm rng tp (Right stuff)+ | tp == "exon" =+ case M.lookup "transcript_id" stuff of+ Just tid+ | tid == c_id xscript -> go gene xscript { c_exons = rng : c_exons xscript } strm++ | otherwise -> Left (GffError fp ln GffIdMismatch) : go gene xscript strm++ Nothing -> Left (GffError fp ln GffParseError) : go gene xscript strm+++ | tp == "transcript" || tp == "cdna" || tp == "mrna" =+ Right xscript :+ case (M.lookup "transcript_id" stuff, M.lookup "gene_id" stuff) of+ (Just tid, Just gid) -> let xscript' = Cdna { c_id = tid+ , c_pos = rng+ , c_gene_id = gid+ , c_gene_symbol = M.lookupDefault "" "gene_name" stuff+ , c_gene_biotype = "" -- XXX+ , c_biotype = "" -- XXX+ , c_description = "" -- XXX+ , c_exons = [] }+ in go gene xscript' strm++ _ -> Left (GffError fp ln GffParseError) : go gene null_cdna strm++ | otherwise = go gene xscript strm+++-- | Parses the random stuff in GFF into a hash table. Returns 'Just+-- (Left _)' if the file uses assignment style ("foo=bar;"), returns+-- 'Just (Right _)' if the file uses statement style ("foo \"bar\";"),+-- otherwise returns Nothing.+parseStuff :: Bytes -> Maybe (Either (M.HashMap Bytes Bytes) (M.HashMap Bytes Bytes))+parseStuff s = Left <$> parse_assignments M.empty s <|>+ Right <$> parse_quoted M.empty s+ where+ parse_assignments !h s0+ | C.null s0 = Just h+ | otherwise = do let (k,s1) = C.break (=='=') s0+ guard . not $ C.null k+ guard . not $ C.null s1+ let (v,s2) = C.break (==';') $ C.tail s1+ parse_assignments (M.insert k v h) (C.drop 1 s2)++ parse_quoted !h s0+ | C.null s0 || C.head s0 == '#' = Just h+ | otherwise = do let (k,s1) = C.break (==' ') s0+ guard . not $ C.null k+ guard $ C.isPrefixOf " \"" s1+ let (v,s2) = C.break (=='"') $ C.drop 2 s1+ guard . not $ C.null s2+ let s3 = C.dropWhile (/=';') s2+ parse_quoted (M.insert k v h) . C.dropWhile isSpace $ C.drop 1 s3+
+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Changelog++## 1.0 (2024-09-03)++- broke twobitreader out of biohazard library
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2024, Udo Stenzel++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Udo Stenzel nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ tool/twobit.hs view
@@ -0,0 +1,122 @@+import Bio.TwoBit+import Bio.TwoBit.Tool+import Codec.Compression.GZip ( decompress )+import Control.Exception+import Control.Monad+import qualified Data.ByteString.Builder as B+import qualified Data.ByteString.Char8 as C+import qualified Data.ByteString.Lazy.Char8 as L+import Options.Applicative+import System.Directory ( renameFile )+import System.IO+import Text.Printf++options :: Parser (IO ())+options = hsubparser $+ command "info" (info tbinfo_options (progDesc "List reference sequences"))+ <> command "tofa" (info tbtofa_options (progDesc "Extract Fasta from 2bit"))+ <> command "fromfa" (info fatotb_options (progDesc "Convert Fasta to 2bit"))+ <> command "fromvcf" (info vcftotb_options (progDesc "Extract reference from dense VCF to 2bit"))+ <> command "transcriptome" (info tbxome_options (progDesc "Transform GTF annotation into a transcriptome"))+++tbinfo_options :: Parser (IO ())+tbinfo_options = mapM_ go <$> some (strArgument (metavar "2BIT-FILE"))+ where+ go f = do ref <- openTwoBit f+ mapM_ (\TBC{..} -> printf "%s\t%d\n" (C.unpack tbc_name) tbc_dna_size) (tbf_chroms ref)+++tbtofa_options :: Parser (IO ())+tbtofa_options = go <$> strArgument (metavar "2BIT-FILE") <*> many (argument rng (metavar "RANGE"))+ where+ rng = do s0 <- str+ if all (/=':') s0+ then pure (s0, Nothing)+ else do (ch, ':':s1) <- pure $ break (':' ==) s0+ (start,'-':s2) : _ <- pure $ reads s1+ if null s2+ then pure (ch, Just (start, Nothing))+ else do (end, "") : _ <- pure $ reads s2+ pure (ch, Just (start,Just end))++ go fp [ ] = do+ ref <- openTwoBit fp+ forM_ (tbf_chroms ref) $ \TBC{..} ->+ putStrLn ('>' : C.unpack tbc_name) >> twoBitToFa maxBound (tbc_fwd_seq 0)++ go fp rns = do+ ref <- openTwoBit fp+ forM_ rns $ \(ch,se) ->+ case (findChrom (C.pack ch) ref, se) of+ (Just tbs, Nothing) -> do printf ">%s\n" ch+ twoBitToFa maxBound $ tbc_fwd_seq tbs 0++ (Just tbs, Just (start,Nothing)) -> do printf ">%s:%d-\n" ch start+ twoBitToFa maxBound $ tbc_fwd_seq tbs start++ (Just tbs, Just (start,Just end)) -> do printf ">%s:%d-%d\n" ch start end+ if start <= end+ then twoBitToFa (end-start) $ tbc_fwd_seq tbs start+ else twoBitToFa (start-end) $ tbc_rev_seq tbs start++ (Nothing, _) -> fail $ "Unknown target: " ++ ch+++tbxome_options :: Parser (IO ())+tbxome_options = go <$> strOption (short 'o' <> long "output" <> metavar "FILE" <> value "-" <> help "Write output to FILE")+ <*> strArgument (metavar "2BIT-FILE")+ <*> many (strArgument (metavar "GFF-File"))+ where+ go :: FilePath -> FilePath -> [FilePath] -> IO ()+ go ofile tbf gffs = do+ ref <- openTwoBit tbf+ withOutputFile ofile $ \h ->+ mapM_ (either (hPutStrLn stderr . displayException) (B.hPutBuilder h . formatCdna ref)) .+ concatMap (uncurry parseAnno)+ =<< readInputs gffs++walkEncodeProgress :: FilePath -> EncodeProgress -> IO ()+walkEncodeProgress fp (Encoded b) = withOutputFile fp $ flip B.hPutBuilder b+walkEncodeProgress fp (EncodeProgress{..}) = do hPrintf stderr "%s: %d kbases (%d Ns), %dkB encoded\n"+ (show ep_seqname) (ep_position `div` 1000) ep_hardmasked (ep_enclength `div` 1024)+ walkEncodeProgress fp ep_tail++gunzip :: L.ByteString -> L.ByteString+gunzip s = if "\x1f\x8b" `L.isPrefixOf` s then decompress s else s++fatotb_options :: Parser (IO ())+fatotb_options = go+ <$> many (strArgument (metavar "FASTA-FILE"))+ <*> strOption (short 'o' <> long "output" <> metavar "FILE" <> value "-" <> help "Write output to FILE")+ where+ go fs fp = walkEncodeProgress fp . faToTwoBit . L.concat . map gunzip . map snd =<< readInputs fs++vcftotb_options :: Parser (IO ())+vcftotb_options = go+ <$> many (strArgument (metavar "VCF-FILE"))+ <*> strOption (short 'o' <> long "output" <> metavar "FILE" <> value "-" <> help "Write output to FILE")+ where+ go fs fp = walkEncodeProgress fp . vcfToTwoBit .+ map L.toStrict . concatMap L.lines . map gunzip . map snd =<<+ readInputs fs+++withOutputFile :: FilePath -> (Handle -> IO a) -> IO a+withOutputFile "-" k = k stdout+withOutputFile f k = bracket (openBinaryFile (f++".#~#") WriteMode) hClose $ \hdl ->+ k hdl >>= \r -> renameFile (f++".#~#") f >> return r++readInputs :: [FilePath] -> IO [( String, L.ByteString )]+readInputs [] = pure . (,) "stdin" <$> L.getContents+readInputs fs = mapM go fs+ where+ go "-" = (,) "stdin" <$> L.getContents+ go f = (,) (show f) <$> L.readFile f+++main :: IO ()+main = id <=< execParser $+ info (options <**> helper)+ (progDesc "Stores genomes compactly in the 2bit format and allows fast extraction of sequences from them." <>+ header "Compact genome storage" <> fullDesc)
+ twobitreader.cabal view
@@ -0,0 +1,54 @@+Cabal-version: 2.0+Name: twobitreader+Version: 1.0+Synopsis: reader for the 2bit file format+Category: Bioinformatics+Description:+ A library and command line tool for working with 2bit files. 2bit is+ a compact file format for genomes introduced by Jim Kent with his BLAT+ suite in the early 2000s. ++Homepage: https://bitbucket.org/ustenzel/twobittool+License: BSD3+License-File: LICENSE+Extra-Source-Files: CHANGELOG.md++Author: Udo Stenzel+Maintainer: u.stenzel@web.de+Copyright: (C) 2024 Udo Stenzel++Build-type: Simple+Tested-with: GHC == 8.10.7, GHC == 9.4.8, GHC == 9.8.2, GHC == 9.10.1++source-repository head+ type: git+ location: https://bitbucket.org/ustenzel/twobitreader.git++Library+ Exposed-modules: Bio.TwoBit Bio.TwoBit.Tool++ Build-depends: base >= 4.8 && < 4.21,+ bytestring >= 0.10.6 && < 0.13,+ mmap == 0.5.*,+ primitive >= 0.6.1 && < 0.10,+ unordered-containers >= 0.2.5.1 && < 0.3++ Default-Language: Haskell2010++ Default-Extensions: BangPatterns, CPP, OverloadedLists, OverloadedStrings, RecordWildCards+++Executable twobit+ Main-is: twobit.hs+ Hs-source-dirs: tool+ Default-Language: Haskell2010++ Build-depends: base >= 4.8 && < 4.21,+ bytestring >= 0.10.6 && < 0.13,+ directory >= 1.2.2 && < 1.4,+ optparse-applicative >= 0.13 && < 0.19,+ twobitreader == 1.0.*,+ zlib >= 0.6 && < 0.8++ Default-Extensions: BangPatterns, OverloadedStrings, RecordWildCards+