GitHub

________                _____
___  __/______ ___________  /______ ________
__  /_  _  __ `/__  ___/_  __/_  _ \__  ___/
_  __/  / /_/ / _(__  ) / /_  /  __/_  /
/_/     \__,_/  /____/  \__/  \___/ /_/

Faster is a little Common Lisp library for working with FASTQ files. Very often in bioinformatics you want to do something like:

  • Take a FASTQ file as input.
  • Loop over each record in the file and do something with it.
  • Write some kind of result (e.g. line of a BED file) for each record, in no particular order.
  • (Optionally) multithread the record processing while still keeping the output unmangled.

Faster is designed for this and is decently fast. Here's a quick example of using it:

(ql:quickload :faster)
(defun gc-content (sequence)
  (/ (count-if (lambda (base)
                 (or (= base (char-code #\G))
                     (= base (char-code #\C))))
               sequence)
     (length sequence)))
(faster:run
  "test/data/example.fastq"
  :report-progress nil
  :worker-threads 4
  :work-function (lambda (record)
                   (cons (faster:parse-sequence-id (faster:id record))
                         (gc-content (faster:seq record))))
  :output-function (lambda (result)
                     (destructuring-bind (id . gc-content) result
                       (format t "Read ~A has GC content ~,3F.~%" id gc-content))))
;; Output:
;;
;; Read example-0 has GC content 0.547.
;; Read example-2 has GC content 0.415.
;; Read example-1 has GC content 0.457.

Basic Interface

(faster:map-fastq function stream) will read successive FASTQ records from stream, parse each into a faster:fastq-record struct, and call function on the struct. stream must be an input stream with element type (unsigned-byte 8).

faster:map-fastq uses an internal buffer to read large chunks of the FASTQ instead of reading line-by-line to try to avoid too many context switches. By default this buffer is rather large (10 megabytes) because our lab works primarily with Oxford Nanopore data, which can have reads that reach megabase length. If you're working with FASTQ files with shorter entries and want to save a bit of RAM you can decrease the buffer size (though really if a single 10mb buffer is a problem in your environment then Common Lisp itself is also likely a problem).

faster:fastq-record is a struct with three slots:

  • faster:id: the sequence identifier.
  • faster:seq: the sequence (i.e. the bases).
  • faster:qs: the quality scores.

All of these slots are of type (simple-array (unsigned-byte 8) (*)), and each slot has an accessor with the same name (e.g. (faster:qs my-fastq-record)).

No extra processing is done on the bytes read from the FASTQ — if you want to turn the sequence into a string or the quality scores into Phred scores you'll need to do that as a separate step. One of the goals of Faster is to make processing FASTQ files fast, and it's often faster to work directly with the bytes and not do any unneeded processing. Faster provides several utility functions for these common transformations:

  • faster:parse-sequence-id
  • faster:bytes-to-quality-scores
  • faster:n-bytes-to-quality-scores

Multithreading Interface

Faster also provides a (faster:run …) function that will use lparallel to distribute reads for parallel processing across worker threads.

(faster:run) takes several arguments:

(filenames &key
 work-function
 output-function
 (interactive t)
 (report-progress t)
 (worker-threads 1)
 (queue-size (* worker-threads 10)))

filenames is a list of file paths (or, for convenience, a bare single path) to process. A path of "-" will read from *standard-input*.

work-function is required and must be a function that takes a single faster:fastq-record, does whatever you want, and returns a result. This function may be invoked in parallel across multiple threads (on different faster:fastq-records).

output-function is required and must be a function that takes the result returned by the work-function and outputs whatever you want. This function will be called on results one-at-a-time, not in parallel, so it's safe to e.g.: write a line to *standard-output* without worrying that multiple reads' output will be mixed together.

interactive determines what happens when an unhandled error occurs in one of the threads. When interactive is t, the normal error-handling of your Lisp is left intact, likely dropping you into a debugger. Passing nil here will cause the entire process to exit with a non-zero exit code (at least on SBCL, behaviour may vary depending on how your implementation treats its exit function when multiple threads are running). This can be useful if you build a binary that is run in the middle of a pipeline — if an error happens you probably want the program to just exit rather than hanging forever.

If report-progress is true, read IDs will be printed to *error-output* as processing progresses. Care is taken to print these sequentially so the output isn't mangled. If nil, no progress is printed.

worker-threads is the number of worker threads to spawn. In addition to the worker threads (i.e. those that call work-function) two or three additional threads will be spawned: an input reader, progress writer (when enabled), and output writer. So e.g. :worker-threads 5 will cause 5 + 1 + 1 + 1 = 8 total threads to be spawned.

queue-size is the size of the input and output queues. Larger queues will allow enqueueing more reads/results to hopefully avoid threads having to stall waiting for input, at the cost of more memory usage.

run will wait until all processing has been completed before it returns, and it does not return any values.

Read the original on github.com ↗