bio-0.5.1: examples/RSelectPE.hs

{- | Randomly select paired end (illumina) reads -}

import System.Environment (getArgs)
import Data.List (isSuffixOf)
import System.Random
import System.IO
import Bio.Sequence

main = do
  (n:as) <- getArgs
  let fs = pairs sufchk as
      p  = read n :: Double
  rs <- map (<p) `fmap` randomRs (0,1.0) `fmap` newStdGen
  s1 <- concat `fmap` mapM readIllumina (map fst fs)
  s2 <- concat `fmap` mapM readIllumina (map snd fs)
  let out = [(s,t) | (s,t,True) <- zip3 s1 s2 rs]
  h1 <- openFile "out.1.txt" WriteMode 
  h2 <- openFile "out.2.txt" WriteMode 
  mapM_ (writer h1 h2) out
  hClose h1
  hClose h2
  
writer :: Handle -> Handle -> (Sequence Nuc,Sequence Nuc) -> IO ()
writer h1 h2 (s,t) | hdchk s t = do
  hWriteIllumina h1 [s]
  hWriteIllumina h2 [t]

pairs :: Show a => (a -> a -> Bool) -> [a] -> [(a,a)]  
pairs check (x:y:zs) 
  | check x y = (x,y):pairs check zs
  | otherwise = error ("Suffix check failed for "++show x++" and "++show y)
pairs _ [] = []
pairs _ _x = error ("Uneven number of args, last was" ++ show _x)

sufchk x y = ".1.txt" `isSuffixOf` x && ".2.txt" `isSuffixOf` y

hdchk x y = f x == f y
  where f = takeWhile (/='/') . toStr . seqlabel