Get rust kmer/fastx code working

Author: Roderick Bovee

Cargo.toml

@@ -0,0 +1,19 @@
+[package]
+name = "needletail"
+version = "0.1.0"
+authors = ["Roderick Bovee <[email protected]>"]
+
+[features]
+default = ["gz"]
+gz = ["flate2"]
+
+[dependencies]
+flate2 = { version="0.2", optional=true }
+memchr = "0.1.11"
+
+[dev-dependencies]
+bencher = "0.1.1"
+
+[[bench]]
+name = "benchmark"
+harness = false

LICENSE

@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2016 Roderick Bovee
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,52 @@
+[![Circle CI](https://circleci.com/gh/onecodex/needletail.svg?style=shield&circle-token=65c2b7d87452dba5e8e3e967133311af478632a4)](https://circleci.com/gh/onecodex/needletail)
+
+# Needletail
+
+Needletail is a MIT-licensed, minimal-copying FASTA/FASTQ parser and k-mer processing library.
+
+The goal is to write a fast *and* well-tested set of functions that more-specialized bioinformatics programs can use.
+Needletail's goal is to be as fast as the `readfq` C library at parsing FASTX files and much (i.e. 25 times) faster than equivalent Python implementations at k-mer counting.
+
+For example, a simple Needletail script can count all the bases in a 2.1 gigabyte FASTQ file in 2.75 seconds while a comparable parser with `readfq` takes 3.52 seconds (see `bench` folder; measured with `%timeit -r 3 -n 3` on an Early 2015 MacBook Pro).
+
+# Example
+
+```rust
+extern crate needletail;
+ues std::env;
+use needletail::{fastx};
+
+fn main() {
+  let filename: String = env::args().nth(1).unwrap();
+
+  let mut n_bases = 0;
+  fastx::fastx_file(&filename[..], |seq| {
+    // seq.0 is the name of the record
+    // seq.1 is the base sequence
+    n_bases += seq.1.len();
+    // seq.2 is an optional quality score
+  });
+  println!("There are {} bases in your file.", n_bases);
+}
+```
+
+# Installation
+
+Needletail requires `rust` and `cargo` to be installed.
+Please use either your local package manager (`homebrew`, `apt-get`, `pacman`, etc) or install these via [rustup](https://www.rustup.rs/).
+
+To install needletail itself for development:
+```shell
+git clone https://github.com/bovee/needletail
+cargo test  # to run tests
+```
+
+# Getting Help
+
+Questions are best directed as GitHub issues.
+
+Hopefully I'll compile the documentation and put it up as a webpage soon too.
+
+# Contributing
+
+Please do! We're happy to discuss/mentor possible additions and/or accept pull requests.

bench/benchmark.c

@@ -0,0 +1,27 @@
+// using https://github.com/lh3/readfq for parsing
+// derived in part from http://stackoverflow.com/questions/19390245/how-to-parse-a-fasta-file-using-kseq-h
+#include <unistd.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "./readfq/kseq.h"
+KSEQ_INIT(FILE*, read)
+
+int main(int argc, char* argv[]) {
+	FILE* fp = fopen(argv[1], "r");
+	if (fp == 0) {
+		perror("fopen");
+		exit(1);
+	}
+	kseq_t *seq = kseq_init(fileno(fp));
+
+	unsigned long slen = 0;
+	while (kseq_read(seq) >= 0) {
+		slen += seq->seq.l;
+	}
+
+	printf("%lu\n", slen);
+	
+	kseq_destroy(seq);
+	fclose(fp);
+	return (0);
+}

bench/benchmark.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+"""
+Unscientific benchmarking of this versus the --release rust
+implementation below using the %timeit Ipython magic (times in sec)
+
+    n_kmers,  py_runtime, rust_runtime
+    6594204,  14.4,       0.578
+
+Both give identical counts on the files tested (and printing kmers out
+and diff'ing the results gives no difference)
+"""
+from __future__ import print_function
+import sys
+
+from Bio.SeqIO import parse
+from Bio.Seq import reverse_complement
+
+
+def slid_win(seq, size=4, overlapping=True):
+    """Returns a sliding window along self.seq."""
+    itr = iter(seq)
+    if overlapping:
+        buf = ''
+        for _ in range(size):
+            buf += next(itr)
+        for l in itr:
+            yield buf
+            buf = buf[1:] + l
+        yield buf
+    else:
+        buf = ''
+        for l in itr:
+            buf += l
+            if len(buf) == size:
+                yield buf
+                buf = ''
+
+filename = sys.argv[1]
+
+n_total = 0
+n_canonical = 0
+
+for s in parse(filename, 'fastq'):
+    uppercase_seq = str(s.upper().seq)
+    for kmer in slid_win(uppercase_seq, 4):
+        canonical = min(kmer, reverse_complement(kmer))
+        if canonical == 'CAGC':
+            n_canonical += 1
+        n_total += 1
+        # print(canonical)
+
+print(n_total, n_canonical)

bench/benchmark.rs

@@ -0,0 +1,20 @@
+extern crate needletail;
+
+use std::env;
+use needletail::{fastx, kmer};
+
+
+fn count_bases(filename: String) {
+    let mut n_bases = 0;
+    fastx::fastx_file(&filename[..], |seq| {
+        n_bases += seq.1.len();
+    }).unwrap();
+    println!("{:?}", n_bases);
+}
+
+
+fn main() {
+    let filename: String = env::args().nth(1).unwrap();
+
+    count_bases(filename);
+}

benches/benchmark.rs

@@ -0,0 +1,53 @@
+#[macro_use]
+extern crate bencher;
+extern crate needletail;
+
+use needletail::{fastx, kmer, bitkmer};
+use bencher::Bencher;
+
+fn check_kmer_speed(bench: &mut Bencher) {
+    let filename = String::from("tests/data/28S.fasta");
+    let ksize = 31;
+
+    bench.iter(|| {
+        let mut n_total = 0;
+        let mut n_canonical = 0;
+        fastx::fastx_file(&filename[..], |seq| {
+            // let normalized_seq = kmer::normalize(Cow::Borrowed(&seq.1), true);
+            for k in seq.1.windows(ksize) {
+                let l = kmer::canonical(k).into_owned();
+                if l == k {
+                    n_canonical += 1;
+                }
+                n_total += 1;
+            }
+        });
+        assert_eq!(213703, n_total);
+        assert_eq!(108521, n_canonical);
+    })
+}
+
+fn check_bitkmer_speed(bench: &mut Bencher) {
+    let filename = String::from("tests/data/28S.fasta");
+    let ksize = 31;
+
+    bench.iter(|| {
+        let mut n_total = 0;
+        let mut n_canonical = 0;
+        fastx::fastx_file(&filename[..], |seq| {
+            for k in bitkmer::BitKmerIter::new(&seq.1, ksize) {
+                let l = bitkmer::canonical(k);
+                if l == k {
+                    n_canonical += 1;
+                }
+                n_total += 1;
+            }
+        });
+        assert_eq!(213703, n_total);
+        assert_eq!(108521, n_canonical);
+    })
+}
+
+
+benchmark_group!(benches, check_kmer_speed, check_bitkmer_speed);
+benchmark_main!(benches);

circle.yml

@@ -0,0 +1,15 @@
+checkout:
+  post:
+    - git config --global --unset url.ssh://[email protected]:.insteadof
+
+dependencies:
+  override:
+    - curl -sS https://static.rust-lang.org/rustup.sh > rustup.sh && chmod +x ./rustup.sh && ./rustup.sh --yes
+
+  cache_directories:
+    - "~/.cargo"
+
+test:
+  override:
+    - cargo build -v --color never
+    - cargo test -v --color never