diff --git a/Bio/TwoBit.hs b/Bio/TwoBit.hs new file mode 100644 --- /dev/null +++ b/Bio/TwoBit.hs @@ -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 + diff --git a/Bio/TwoBit/Tool.hs b/Bio/TwoBit/Tool.hs new file mode 100644 --- /dev/null +++ b/Bio/TwoBit/Tool.hs @@ -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 + diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# Changelog + +## 1.0 (2024-09-03) + +- broke twobitreader out of biohazard library diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/tool/twobit.hs b/tool/twobit.hs new file mode 100644 --- /dev/null +++ b/tool/twobit.hs @@ -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) diff --git a/twobitreader.cabal b/twobitreader.cabal new file mode 100644 --- /dev/null +++ b/twobitreader.cabal @@ -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 +