matsen/pplacer — Issue #233: guppy nbc
Author: matsen
Created: 2012-02-28 14:36:12+00:00
State: open
Labels: top

nbc = Naive Bayes Classifier.

The code here should be factored as usual, because eventually we are going to be doing an NBC + phylogenetic classification hybrid.

Note that there is some chance in the future we will want to change to a 32-bit int for `freq_table`, or even a float. I will use "vector" to mean Bigarray.Array1. There is some reasonable code to play with these things in Gsl_vector.

Flags include:

-c Ref package
--word-length The length of the word used to do the NBC classification (Default 8)
--target-rank The desired most specific rank for classification (Default genus)

Let `n_words` = 4^`word_length`. Assert that `n_words < max_int`.

We will denote the corresponding variables with underscores in place of dashes.

Start by grouping the reference sequences at the given target rank, discarding those that do not have a classification at that rank. Say there are `n_taxids` taxids represented by at least one sequence at the target rank.

Allocate an int bigarray.Array2 `freq_table` of width `n_words` and height `n_taxids` with all zero entries.

We will need a subroutine `word_to_int` to turn an `[ACGT]{word_length}` into an integer. This will work by turning every A, C, G, and T into a pair of bits like

A -> 0 0
C -> 0 1
G -> 1 0
T -> 1 1

Given a word, we concatenate these bits together into an integer. Please make an inverse function and then test round-tripping.

Using this routine, build another `add_counts_by_seq` that takes a sequence and an int big-vector of length `n_words` and marches along its de-gapped version, finding words and adding to the corresponding entry. Specifically, for every word W of the ungapped sequence, increase the corresponding entry of the vector by one. For the time being, ignore words that are not exclusively composed of {A,C,G,T}.

Using this routine, fill out the entries in the `freq_table` as follows. If a given sequence has taxonomic classification X, pull out the row for that taxid using slice and pass it along the sequence to `add_counts_by_seq`.

Make an float vector `prior_counts` such that the entry for word i is (n(w_i) + 0.5) / (N + 1), where n(w_i) is the total number of observations of w_i and N is the total number of sequences.

Then build a float Array2 called `taxid_word_counts` that is `n_taxids` by `n_words`, where the entry at (x,i) is log (m(w_i) + `prior_counts`[i]) - `denom`; m(w_i) is the number of sequences with word i for taxid x and `denom` = log (M + 1) where M is the number of sequences with taxid x.

Then, to classify a sequence, use `add_counts_by_seq` to get a vector of k-mer counts for that sequence. Multiply this vector by the `taxid_word_counts` matrix using `Linear_utils.mat_vec_mul` (into a pre-allocated vector) and find the index of the maximum value using `Gsl_vector.max_index`. The taxid corresponding to this index is our classification.

This classification should then get output via the previous classification machinery.