kmer-boolean

This utility tests if a specified kmer is or is not in a set of FASTA sequences provided on standard input, for a given k, returning the according "true" or "false" result. Alternatively, the binary can instead report kmers that are found (--present) or not found (--absent), or either case (--all).

Memory usage

Internally, this test keeps an array of bits to minimize the memory overhead of storing per-kmer presence or absence state. This requires at least 2^2k-3 bytes to store the bitarray. Querying 16mers, for example, will require 537 MB of memory.

For the C++ binary, an additional 8 MB buffer is reserved for storing intermediate sequence data that streams in from the input FASTA file. If the --read-in-all-at-once option is used, the sequence data is read into memory all at once. It is strongly recommended to instead use the default streaming option to minimize memory usage and runtime overhead required to read all sequences into memory.

Runtime

To explore runtime characteristics, we used hg38 assembly data, looking for kmers that are absent in this genome build from k=2 up to 12, to get a general sense of how runtimes trend.

This test is naturally parallelized by chromosome, as any kmer absent in one chromosome is absent over all chromosomes, by definition. We include scripts compatible with a Slurm job scheduler for local testing.

While we include scripts to parallelize the generation of runtime data for all k-by-chromosome combinations, plotting uses the maximum runtime among all chromosomes as the natural bound for querying kmers for the given value of k.

Some work is still left to do to optimize data structures used for bitsetting kmers.

Notes regarding FASTA input

Note that searches are not performed on "canonical" DNA kmers, but on all unique kmers. For example, a search for the AG 2mer will be treated separately from the reverse complement CT 2mer.

The input FASTA file may contain one or more records (so-called "multi-FASTA"). Each FASTA record in a multi-FASTA file is scanned separately, but each record contributes to the overall kmer query report. Split a multi-FASTA file if you want to query kmer distributions for individual records.

Kmers containing hard-masked bases (N) are ignored. Soft-masked bases (lowercase bases) are included in queries.

An 8 MB buffer is kept of sequence data streaming in from the input FASTA. It is also possible to use the --read-in-all-at-once option to read the entire FASTA records into memory. Using this option would not be recommended for genome-scale input files.

C++

This C++ implementation includes a custom bitset container which can be sized at runtime. The STL bitset library can only be sized at compile time and thus is not used here.

Compilation

Run make to build the kmer-boolean binary. Run make clean to clean up temporary files.

Usage

Options

Run kmer-boolean --help to get a listing of options:

kmer-boolean
  version: 1.1
  author:  Alex Reynolds

  Usage:

  $ kmer-boolean [arguments] < input

  Test if specified kmer is or is not in a set of
  FASTA sequences provided on standard input, for a
  given k, returning the according 'true' or 'false'
  result.

  General Options:

  --k=n                          K-value for kmer length (integer)
  --query-kmer=s                 Test or query kmer (string)
 [--present | --absent | --all]  If query kmer is omitted, print all
                                 present or absent kmers, or all kmers
  --read-in-all-at-once          Read all sequence data into memory
                                 before processing (not recommended for
                                 very large FASTA inputs; default is to
                                 stream through input in chunks)
                                 
  Process Flags:

  --help                         Show this usage message
  --version                      Show binary version

Examples

You can return a list of all kmers found in the input via --present:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-boolean --k=2 --present
AC found
AT found
AG found
CA found
CC found
CT found
CG found
TA found
TC found
TT found
GA found
GG found

Conversely, you can return a list of kmers not found in the input via --absent:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-boolean --k=2 --absent
AA not found
TG not found
GC not found
GT not found

You can also use --all to return a list of all kmers, whether found or not found in the input FASTA:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-boolean --k=3 --all
AAA not found
AAC not found
AAT not found
AAG not found
ACA not found
ACC found
ACT not found
...
GTT not found
GTG not found
GGA found
GGC not found
GGT not found
GGG found

If a filter option (--present, etc.) is not specified, kmer-boolean will return a list of all present kmers.

You may also query for a specific kmer with the --query-kmer=<kmer> option:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-boolean --k=5 --query-kmer="ACGGC"
ACGGC not found

If the queried kmer is found, the exit code for the application is zero ("success"). If the kmer is not found, the exit code is non-zero ("failure").

Python

The kmer-boolean.py script uses the Python bitarray library to set up the required number of booleans for the specified k, and tests the boolean value for the specified kmer, or else returns those kmers which are observed in the input.

Note that searches are performed not on "canonical" kmers, but all unique kmers. In other words, for example, a search for the AG 2mer will be treated separately from the reverse complement CT 2mer.

Testing

This has been tested under Python 3.6.3 on OS X. Use of conda or pip may be required for installation of the bitstring library. The makefile offers a python_test target to run expected "true" and "false" tests of searching for a pair of 2mers from some ad-hoc sequences piped in via standard input to the kmer-boolean.py script.