pipes-illumina (empty) → 0.1.0.0
raw patch · 5 files changed
+228/−0 lines, 5 filesdep +basedep +bytestringdep +directorysetup-changed
Dependencies added: base, bytestring, directory, filepath, pipes, pipes-bgzf
Files
- LICENSE +30/−0
- Pipes/Illumina.hs +126/−0
- Pipes/NextSeq.hs +48/−0
- Setup.hs +2/−0
- pipes-illumina.cabal +22/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014-2015, Rob O'Callahan++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 Your name here 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.
+ Pipes/Illumina.hs view
@@ -0,0 +1,126 @@+{-# LANGUAGE BangPatterns #-}+module Pipes.Illumina where++import Data.Bits+import Data.Word+import Data.Int+import Foreign.Ptr+import Foreign.Storable+import Foreign.ForeignPtr.Unsafe+import Foreign.ForeignPtr.Safe+import Foreign.Marshal.Alloc+import Data.ByteString.Internal+import Control.Applicative+import Control.Monad+import System.IO+import qualified Data.ByteString as B+import Data.ByteString.Unsafe+import Pipes+import Pipes.Bgzf++bclBgzfProducer :: MonadIO m => [Handle] -> Producer (ByteString, ByteString) m ()+bclBgzfProducer [] = return ()+bclBgzfProducer hdls = start where+ start = do+ mblocks <- getBlocks+ case mblocks of+ Nothing -> return ()+ Just bs -> do+ b2fPipe $ map (B.drop 4) bs+ go+ go = do+ mblocks <- getBlocks+ case mblocks of+ Nothing -> return ()+ Just bs -> do+ b2fPipe bs+ go+ getBlocks = do+ mbchunks <- forM hdls $ \h -> do+ header <- liftIO $ B.hGet h 18+ case parseHeader header of+ 28 -> return Nothing+ blocklen -> do+ chunk <- liftIO $ B.hGet h (blocklen - 18)+ return $ Just chunk+ return $ map inflateBlock <$> sequence mbchunks+ cycs = length hdls+ b2fPipe chunks = do+ let l = B.length $ head chunks+ forM_ [0..l-1] $ \i -> do+ r <- liftIO $ do+ seqfp <- mallocByteString cycs+ qualfp <- mallocByteString cycs+ let seqp :: Ptr Word8+ seqp = unsafeForeignPtrToPtr seqfp+ qualp = unsafeForeignPtrToPtr qualfp+ forM_ (zip chunks [0..]) $ \(c,j) -> do+ case unsafeIndex c i of+ 0 -> do+ pokeElemOff seqp j 78 + pokeElemOff qualp j 35+ w -> do+ pokeElemOff seqp j $ case w .&. 3 of+ 0 -> 65+ 1 -> 67+ 2 -> 71+ 3 -> 84+ pokeElemOff qualp j $ 33 + w `shiftR` 2+ touchForeignPtr seqfp+ touchForeignPtr qualfp+ let !sq = PS seqfp 0 cycs+ !qual = PS qualfp 0 cycs+ return $! (sq, qual)+ yield r++filterProducer :: Handle -> Producer Bool IO ()+filterProducer hdl = go where+ numclusters :: IO Int32+ numclusters = do+ buf <- mallocBytes 12+ 12 <- hGetBufSome hdl buf 12+ cnt <- peek (buf `plusPtr` 8)+ free buf+ return cnt+ go = do+ clusts <- liftIO $ fromIntegral <$> numclusters+ buf <- liftIO $ mallocBytes 1+ replicateM_ clusts $ do+ 1 <- liftIO $ hGetBufSome hdl buf 1+ b <- liftIO $ peek buf+ yield b+ liftIO $ free buf++locsProducer :: Handle -> Producer (Float, Float) IO ()+locsProducer hdl = go where+ numclusters :: IO Int32+ numclusters = do+ buf <- mallocBytes 12+ 12 <- hGetBufSome hdl buf 12+ cnt <- peek (buf `plusPtr` 8)+ free buf+ return cnt+ go = do+ clusts <- liftIO $ fromIntegral <$> numclusters+ buf <- liftIO $ mallocBytes 8+ replicateM_ clusts $ do+ 8 <- liftIO $ hGetBuf hdl buf 8+ x <- liftIO $ peek buf+ y <- liftIO $ peekElemOff buf 1+ yield (x, y)+ liftIO $ free buf++bciProducer :: Handle -> Producer (Int, Int) IO ()+bciProducer hdl = go where+ go = do+ buf <- liftIO $ (mallocBytes 8 :: IO (Ptr Int32))+ let go' = do+ r <- liftIO $ hGetBuf hdl buf 8+ case r of+ 8 -> do+ a <- liftIO $ peek buf+ b <- liftIO $ peekElemOff buf 1+ yield (fromIntegral a, fromIntegral b)+ go'+ 0 -> liftIO $ free buf+ go'
+ Pipes/NextSeq.hs view
@@ -0,0 +1,48 @@+{-# LANGUAGE OverloadedStrings #-}++module Pipes.NextSeq where++import Control.Applicative+import Control.Monad+import System.IO+import qualified Data.ByteString as B+import Data.ByteString (ByteString)+import Data.ByteString.Unsafe+import Data.List+import Pipes+import qualified Pipes.Prelude as P+import Pipes.Illumina+import System.FilePath+import System.Directory++data Cluster = Cluster+ { sq :: !ByteString+ , qual :: !ByteString+ , xCoord :: !Float+ , yCoord :: !Float+ , passedFilter :: !Bool+ , lane :: !Int+ , tile :: !Int }++nextSeqProducer :: FilePath -> Producer Cluster IO ()+nextSeqProducer rundir = forM_ [1..4] proclane where+ proclane ln = do+ let lanestr = "L00" ++ show ln+ bcldir = joinPath [rundir, "Data", "Intensities", "BaseCalls", lanestr]+ allconts <- liftIO $ getDirectoryContents bcldir+ let bcls = sort $ filter (\q -> takeExtension q == ".bgzf") allconts+ locshdl <- liftIO $ openFile (joinPath [rundir, "Data", "Intensities", lanestr, "s_" ++ show ln ++ ".locs"]) ReadMode+ filthdl <- liftIO $ openFile (joinPath [bcldir, "s_" ++ show ln ++ ".filter"]) ReadMode+ bcihdl <- liftIO $ openFile (joinPath [bcldir, "s_" ++ show ln ++ ".bci"]) ReadMode+ bclhdls <- liftIO $ mapM (\r -> openFile (joinPath [bcldir, r]) ReadMode) bcls+ let bclprod = bclBgzfProducer bclhdls+ tileprod = for (bciProducer bcihdl) $ \(tnum, cnt) -> replicateM_ cnt $ yield tnum+ filtprod = filterProducer filthdl+ locsprod = locsProducer locshdl+ zipd = P.zip tileprod $ P.zip filtprod $ P.zip locsprod bclprod+ for zipd $ \(tnum, (pf, ((x,y), (sq, qu)))) ->+ yield $ Cluster sq qu x y pf ln tnum+ liftIO $ hClose locshdl+ liftIO $ hClose filthdl+ liftIO $ hClose bcihdl+ liftIO $ mapM_ hClose bclhdls
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ pipes-illumina.cabal view
@@ -0,0 +1,22 @@+name: pipes-illumina+version: 0.1.0.0+synopsis: Illumina NGS data processing+description: Streaming of Bcl and other Illumina file formats+homepage: http://github.com/rcallahan/pipes-illumina+license: BSD3+license-file: LICENSE+author: Rob O'Callahan+maintainer: ropoctl@gmail.com+copyright: 2014-2015+category: Codec+build-type: Simple+cabal-version: >=1.10++library+ exposed-modules: Pipes.Illumina, Pipes.NextSeq+ build-depends: base >= 4.7 && < 5, bytestring, directory, filepath, pipes, pipes-bgzf > 0.2+ default-language: Haskell2010++source-repository head+ type: git+ location: https://github.com/rcallahan/pipes-illumina