aboutsummaryrefslogtreecommitdiff
path: root/xenomapper.hs
blob: 4546ab2bc364abaee75152220dd3d0f1d7cd8de7 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ViewPatterns #-}
module Main where

import Bio.Prelude hiding ((.))
import Bio.Bam
import Control.Monad
import Control.Monad.Log
import Control.Monad.IO.Class
import System.Environment
import qualified Streaming.Prelude as S
import Streaming.Prelude (Stream, Of)

filterBams :: (MonadIO m, MonadLog m, MonadMask m) => Stream (Of BamRaw) m () ->  Stream (Of BamRaw) m () ->  Stream (Of BamRaw) m ()
filterBams a b = do
  str0 <- lift $ S.next a
  str1 <- lift $ S.next b
  loop str0 str1
  where
    loop a@(Right (a', a'')) b@(Right (b', b'')) =
      case compare (qname a') (qname b') of
        EQ -> do
          when (as a' > as b') $ S.yield a'
          filterBams a'' b''
        LT -> do
          n <- lift $ S.next a''
          loop n b
        GT -> do
          n <- lift $ S.next b''
          loop a n
          
    loop _ _ = return ()
    
    qname = b_qname . unpackBam
    as = extAsInt 0 "AS" . unpackBam

main = do
  args <- getArgs
  if length args /= 3 then do
    path <- getProgName
    hPutStrLn stderr $ unwords [ "usage:", path, "output.bam", "input1.bam", "input2.bam" ]
  else
    withLogging_ (LoggingConf Warning Error Error 10 True) $ do
      decodeBamFiles (tail args) (\[(hdr, prim), (_, sec)] -> writeBamFile (head args) hdr $ filterBams prim sec)