1
\$\begingroup\$

I'm a bioinformatician and mainly use python, but I'm trying to learn Haskell. So far I like this language because it's beautiful and fun, but my Haskell scripts usually end up being much slower than my Python equivalent.

Here is my latest attempt:

-- Define type synonyms for fields in a fastq record.
type Name = String
type Nucleotides = String
type Qualities = String

data FastqRecord = FastqRecord Name Nucleotides Qualities
{-
Show is implemented so that the string representation of a FastqRecord
is actually its fastq formatted representation.
-}
instance Show FastqRecord where
    show (FastqRecord name s q) = unlines [name, s, "+", q]

-- Parse lines four by four, ignore the third of each batch.
getFastqRecords :: [String] -> [FastqRecord]
getFastqRecords [] = []
getFastqRecords (l1 : l2 : _ : l4 : ls) = FastqRecord l1 l2 l4 : getFastqRecords ls

{-
We assume that the records come in sorted by sequence.
Successive records having the same sequence are fused,
arbitrarilly keeping the first name, and taking the
highest quality at each position.
-}
fuseFastqRecords :: [FastqRecord] -> [FastqRecord]
fuseFastqRecords [] = []
fuseFastqRecords [r] = [r]
fuseFastqRecords (r1@(FastqRecord n s1 q1) : r2@(FastqRecord _ s2 q2) : rs)
    | s1 == s2 = fuseFastqRecords (FastqRecord n s1 (bestquals q1 q2) : rs)
    | otherwise = r1 : fuseFastqRecords (r2 : rs)

{-
For each position, we want to keep the highest quality
among the records that had the same sequence.
In fastq format, qualities are integers represented as ascii characters
Using max on Char does what we want.
-}
bestquals :: Qualities -> Qualities -> Qualities
bestquals = zipWith max

main :: IO ()
main = interact (concatMap show . fuseFastqRecords . getFastqRecords . lines)

On very small test input, this scripts runs 10 times faster than my python implementation.

On "real" input, the python code (using pypy) uses much more ram than the Haskell version, but runs in half the time.

If I understood my readings so far, Strings in Haskell are not implemented efficiently, so I guess I would have to use another type to process my text. But the representation as a list of Char is actually very convenient for one of the operations I want to do (see the comments about the bestquals function).

What would you recommend to do?

\$\endgroup\$

2 Answers 2

1
\$\begingroup\$

I found out that using Data.ByteString.Lazy.Char8 could significantly improve the running time (about one third of the one I posted initially).

Here is the code, taking into account some suggestions by Gurkenglas:

import qualified Data.ByteString.Lazy.Char8 as C

-- Define type synonyms for fields in a fastq record.
type Name = C.ByteString
type Nucleotides = C.ByteString
type Qualities = C.ByteString

data Fastq = Fastq Name Nucleotides Qualities

formatFastq :: Fastq -> C.ByteString
formatFastq (Fastq n s q) = C.unlines [n, s, C.pack "+", q]

-- Parse lines four by four, ignore the third of each batch.
getFastqs :: [C.ByteString] -> [Fastq]
getFastqs [] = []
getFastqs (l1 : l2 : _ : l4 : ls) = Fastq l1 l2 l4 : getFastqs ls

{-
We assume that the records come in sorted by sequence.
Successive records having the same sequence are fused,
arbitrarilly keeping the first name, and taking the
highest quality at each position.
-}
fuseFastqs :: [Fastq] -> [Fastq]
fuseFastqs [] = []
fuseFastqs [r] = [r]
fuseFastqs (r1@(Fastq n s1 q1) : r2@(Fastq _ s2 q2) : rs)
    | s1 == s2 = fuseFastqs (Fastq n s1 (bestquals q1 q2) : rs)
    | otherwise = r1 : fuseFastqs (r2 : rs)

{-
For each position, we want to keep the highest quality
among the records that had the same sequence.
In fastq format, qualities are integers represented as ascii characters
Using max on Char does what we want, and seems to work on Bytestring elements.
-}
bestquals :: Qualities -> Qualities -> Qualities
bestquals q1 q2 = C.pack (C.zipWith max q1 q2)

processLines :: [C.ByteString] -> [C.ByteString]
processLines ls = map formatFastq (fuseFastqs (getFastqs ls))

main :: IO ()
main = C.interact (C.concat . processLines . C.lines)

It would be nice to be able to avoid the ugly C.pack in bestquals.

\$\endgroup\$
0
\$\begingroup\$

Have you tried compiling with the -O2 flag, as in ghc -O2 yourfile.hs? Other than that, perhaps refactoring the code to use tried and true library functions might help:

import Data.List
import Data.List.Extra -- You'll need to install the extra package.

data Fastq = Fastq
  { name :: String
  , nucleotides :: String
  , qualities :: String
  }

-- Not using type aliases for the Strings is just a matter of taste and doesn't matter to the resulting code.

-- Also by taste, no need to say it's a record because that's what the data definition does. Unless that term also exists in bioinformatics?

fuseFastqs :: [Fastq] -> [Fastq]
fuseFastqs = map fuse . groupOn nucleotides where
  fuse :: [Fastq] -> Fastq
  fuse xs@((Fastq n s _):_) = Fastq n s $ map maximum $ transpose $ map nucleotides xs

By the way, the closest we've got to a typeclass law for Show is that it produces strings that you can paste into Haskell code to recover the shown value. Pretty-printing or format conversion functions are usually called something else.

\$\endgroup\$
1
  • \$\begingroup\$ The -O2 flag seemed to very slightly reduce the running time. I tried to use ByteString instead of String, and this reduced the running time to one third, that is, better than the Python version. I then tried to adapt your version of fuseFastqs, but on a first attempt, the program ran out of memory. On a second attempt, the running time was about the same as my python version. It may well be that my way of combining your approach with a ByteString converted approach is far from being optimally done. Thanks for the "taste" remarks. \$\endgroup\$
    – bli
    Commented Aug 2, 2016 at 15:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.