2.6 Machine learning in bioinformatics ¶

2.6.3 scikit-learn¶

In this chapter we'll implement several machine learning classifiers so we can gain an in-depth understanding of how they work. In practice though, there are many mature machine learning libraries that you'd want to use. scikit-learn is a popular and well-documented Python library for machine learning which many bioinformatics researchers and software developers use in their work.

2.6.4 Defining the problem¶

We'll explore machine learning classifiers in the context of a familiar topic: taxonomic classification of 16S rRNA sequences. We previously explored this problem in Sequence Homology Searching, so it is likely worth spending a few minutes skimming that chapter if it's not fresh in your mind.

Briefly, the problem that we are going to address here is as follows. We have a query sequence ($q_i$) which is not taxonomically annotated (meaning we don't know the taxonomy of the organism whose genome it is found in), and a reference database ($R$) of taxonomically annotated sequences ($r_1, r_2, r_3, ... r_n$). We want to infer a taxonomic annotation for $q_i$. We'll again work with the Greengenes database, which we'll access using QIIME default reference project. Greengenes is a database of 16S rRNA gene sequences. (This should all sound very familiar - if not, I again suggest that you review Sequence Homology Searching.)

This time, instead of using sequence alignment to identify the most likely taxonomic origin of a sequence, we'll train classifiers by building kmer-based models of the 16S sequences of taxa in our reference database. We'll then run our query sequences through those models to identify the most likely taxonomic origin of each query sequence. Since we know the taxonomic origin of our query sequences in this case, we can evaluate the accuracy of our classifiers by seeing how often they return the known taxonomy assignment. If our training and testing approaches are well-designed, the performance on our tests will inform us of how accurate we can expect our classifier to be on data where the actual taxonomic origin is unknown.

Let's jump in...

2.6.5 Naive Bayes classifiers¶

The first classifier we'll explore is the popular and relatively simple Naive Bayes classifier. This classifier uses Bayes Theorem to determine the most likely label for an unknown input based on a probabilistic model it has constructed from training data. (The preceding text needs work.) The model that is constructed is based on user-defined features of the sequences. The most commonly used features for sequence classification tasks such as this is overlapping kmers.

We'll begin by importing some libraries that we'll use in this chapter, and then preparing our reference database and query sequences as we did previously.

In [1]:
%pylab inline

from IPython.core import page
page.page = print

import pandas as pd
import skbio
import numpy as np
import itertools
import collections

Populating the interactive namespace from numpy and matplotlib

In [2]:
import qiime_default_reference as qdr

reference_taxonomy = {}
for e in open(qdr.get_reference_taxonomy()):
seq_id, seq_tax = e.strip().split('\t')
reference_taxonomy[seq_id] = seq_tax

# Load the reference sequences, and associate the taxonomic annotation with
reference_db = []
for e in skbio.io.read(qdr.get_reference_sequences(), format='fasta', constructor=skbio.DNA):
if e.has_degenerates():
# For the purpose of this lesson, we're going to ignore sequences that contain
# degenerate characters (i.e., characters other than A, C, G, or T)
continue
reference_db.append(e)

print("%s sequences were loaded from the reference database." % len(reference_db))

88452 sequences were loaded from the reference database.

In [3]:
reference_db[0]

Out[3]:
DNA
-----------------------------------------------------------------------
'description': ''
'id': '1111883'
'taxonomy': 'k__Bacteria; p__Gemmatimonadetes; c__Gemm-1; o__; f__;
g__; s__'
Stats:
length: 1428
has gaps: False
has degenerates: False
has non-degenerates: True
GC-content: 61.90%
-----------------------------------------------------------------------
0    GCTGGCGGCG TGCCTAACAC ATGTAAGTCG AACGGGACTG GGGGCAACTC CAGTTCAGTG
60   GCAGACGGGT GCGTAACACG TGAGCAACTT GTCCGACGGC GGGGGATAGC CGGCCCAACG
...
1320 GCCGCGGTGA ATACGTTCCC GGGCCTTGTA CACACCGCCC GTCACGCCAT GGAAGCCGGA
1380 GGGACCCGAA ACCGGTGGGC CAACCGCAAG GGGGCAGCCG TCTAAGGT
In [4]:
reference_db[-1]

Out[4]:
DNA
----------------------------------------------------------------------
'description': ''
'id': '4483258'
'taxonomy': 'k__Archaea; p__Crenarchaeota; c__Thermoprotei;
o__Thermoproteales; f__Thermoproteaceae; g__; s__'
Stats:
length: 2123
has gaps: False
has degenerates: False
has non-degenerates: True
GC-content: 58.36%
----------------------------------------------------------------------
0    CTGGTTGATC CTGCCGGACC CGACCGCTAT CGGGGTGGGG CTTAGCCATG CGAGTCAAGC
60   GCCCCAGGGA CCCGCTGGGG TGCGGCGCAC GGCTCAGTAA CACGTGGCCA ACCTACCCTC
...
2040 ATAATCTCCT TATTGTCTGA TCCTTATGCA TTTTCCTTTG GCCCATCCCG TGAATACGCG
2100 CGGTGAATAC GTCCCTGCCC CTT

We'll select a random subset of the reference database to work with here.

In [5]:
reference_db = np.random.choice(reference_db, 500, replace=False)
print("%s sequences are present in the subsampled database." % len(reference_db))

500 sequences are present in the subsampled database.


The first thing our Naive Bayes classifier will need is the set of all possible words of length k. This will be dependent on the value of k and the characters in our alphabet (i.e., the characters that we should expect to find in the reference database). This set is referred to as W, and can be computed as follows. Given the following alphabet, how many kmers of length 2 are there (i.e., 2-mers)? How many 7-mers are there? How many 7-mers are there if there are twenty characters in our alphabet (as would be the case if we were working with protein sequences instead of DNA sequences)?

In [6]:
alphabet = skbio.DNA.nondegenerate_chars
k = 2

def compute_W(alphabet, k):
return set(map(''.join, itertools.product(alphabet, repeat=k)))

W = compute_W(alphabet, k)
print('Alphabet contains the characters: %s' % ', '.join(alphabet))
print('For an alphabet size of %d, W contains %d length-%d kmers.' % (len(alphabet), len(W), k))

Alphabet contains the characters: T, A, G, C
For an alphabet size of 4, W contains 16 length-2 kmers.


scikit-bio provides methods for identifying all kmers in a skbio.DNA sequence object, and for computing the kmer frequencies. This information can be obtained for one of our reference sequences as follows:

In [7]:
kmers = reference_db[0].iter_kmers(k=k)
for kmer in kmers:
print(kmer, end=' ')

AA AC CG GA AA AC CG GC CT TG GG GC CG GG GC CG GT TG GC CC CT TA AA AC CA AC CA AT TG GC CA AA AG GT TC CG GA AA AC CG GA AG GA AA AA AG GT TT TT TC CC CT TT TC CG GG GG GG GA AA AT TG GA AG GT TA AA AA AG GT TG GG GC CA AG GA AC CG GG GG GT TG GA AG GT TA AA AC CA AC CG GT TG GG GG GC CA AA AC CC CT TG GC CC CC CG GG GA AG GG GT TG GG GG GG GG GA AT TA AA AC CC CC CT TC CC CG GA AA AA AG GG GA AG GG GG GC CT TA AA AT TA AC CC CG GC CA AT TG GA AC CA AC CC CT TC CG GA AG GC CA AC CT TT TA AG GG GG GG GC CT TC CG GA AG GA AT TC CA AA AA AG GC CC CA AG GC CC CT TC CT TG GC CT TT TG GC CA AC CG GC CT TG GG GT TG GC CC CA AC CC CG GG GA AG GG GG GA AC CC CC CG GC CG GC CC CC CC CA AT TT TA AG GC CT TT TG GT TT TG GG GT TG GA AG GG GT TA AA AC CG GG GC CT TC CA AC CC CA AA AG GG GC CC CG GC CG GA AT TG GG GG GT TA AG GC CC CG GG GC CC CT TG GA AG GA AG GG GG GT TG GT TT TC CG GG GC CC CA AC CA AC CT TG GG GA AA AC CT TG GA AG GA AC CA AC CG GG GT TC CC CA AG GA AC CT TC CC CT TA AC CG GG GG GA AG GG GC CA AG GC CA AG GT TG GA AG GG GA AA AT TT TT TT TG GC CG GC CA AA AT TG GG GG GC CG GA AA AA AG GC CC CT TG GA AC CG GC CA AG GC CG GA AC CG GC CC CG GC CG GT TG GG GA AG GG GA AC CG GA AA AG GG GT TC CT TT TC CG GG GA AT TC CG GT TA AA AA AC CT TC CC CT TG GT TC CA AA AG GC CG GG GG GA AC CG GA AA AT TC CC CT TG GC CG GG GG GG GA AT TG GA AA AT TA AA AG GC CT TC CC CG GC CA AG GT TT TG GA AC CG GG GT TA AC CC CG GC CT TG GG GA AG GG GA AA AG GC CC CC CC CG GG GC CT TA AA AC CT TC CC CG GT TG GC CC CA AG GC CA AG GC CC CG GC CG GG GT TA AA AT TA AC CG GG GA AG GG GG GG GG GC CC CA AG GC CG GT TT TG GT TT TC CG GG GA AA AT TT TA AT TT TG GG GG GC CG GT TA AA AA AG GG GG GC CG GC CG GT TA AG GG GC CG GG GC CC CG GA AG GC CA AA AG GT TC CT TG GA AG GG GT TG GA AA AA AT TC CC CC CC CC CG GG GC CT TC CA AA AC CC CG GG GG GG GA AA AC CT TG GC CC CT TC CG GG GA AA AA AC CT TG GC CT TC CG GG GC CT TT TG GA AG GT TC CC CC CG GG GA AG GA AG GG GG GT TA AG GC CG GG GA AA AT TT TC CC CC CA AG GT TG GT TA AA AC CG GG GT TG GA AA AA AT TG GC CG GT TA AG GA AT TA AT TT TG GG GG GA AG GG GA AA AC CA AC CC CA AG GT TG GG GC CG GA AA AG GG GC CG GG GC CT TA AC CC CT TG GG GA AC CG GG GT TA AC CT TG GA AC CG GC CT TG GA AG GG GC CG GC CG GA AA AA AG GC CT TA AG GG GG GG GA AG GC CA AA AA AC CA AG GG GA AT TT TA AG GA AT TA AC CC CC CT TG GG GT TA AG GT TC CC CT TA AG GC CT TG GT TA AA AA AC CG GA AT TG GG GG GC CA AC CT TT TG GG GT TG GT TT TG GC CG GG GG GT TA AT TC CG GA AC CC CC CC CT TG GC CA AG GT TG GC CC CG GA AA AG GC CT TA AA AC CG GC CA AT TT TA AA AG GT TG GC CC CC CC CG GC CC CT TG GG GG GG GA AG GT TA AC CG GG GT TC CG GC CA AA AG GG GC CT TG GA AA AA AC CT TC CA AA AA AG GG GA AA AT TT TG GA AC CG GG GG GG GG GC CC CC CG GC CA AC CA AA AG GC CG GG GT TG GG GA AG GC CA AT TG GT TG GG GT TT TT TA AA AT TT TC CG GA AC CG GC CA AA AC CG GC CG GA AA AG GA AA AC CC CT TT TA AC CC CT TG GG GG GT TT TT TG GA AA AC CT TG GC CG GG GC CG GG GA AC CA AG GC CT TC CC CA AG GA AG GA AT TG GG GA AG GC CG GT TC CC CC CT TT TC CG GG GG GG GC CT TG GC CC CG GT TG GG GA AG GG GT TG GC CT TG GC CA AT TG GG GC CT TG GT TC CG GT TC CA AG GC CT TC CG GT TG GT TC CG GT TG GA AG GA AT TG GT TT TG GG GG GT TT TA AA AG GT TC CC CC CG GC CA AA AC CG GA AG GC CG GC CA AA AC CC CC CC CT TG GC CC CC CT TT TA AG GT TT TG GC CC CA AG GC CG GG GT TT TC CG GG GC CC CG GG GG GA AA AC CT TC CT TA AG GG GG GG GG GA AC CT TG GC CC CG GG GT TG GA AT TA AA AA AC CC CG GG GA AG GG GA AA AG GG GT TG GG GG GG GA AT TG GA AC CG GT TC CA AA AG GT TC CA AT TC CA AT TG GG GC CC CT TT TT TA AT TG GT TC CC CA AG GG GG GC CT TA AC CA AC CA AC CG GT TG GC CT TA AC CA AA AT TG GG GA AC CG GG GT TA AC CA AA AA AG GG GG GC CT TG GC CG GA AT TA AC CC CG GC CG GA AG GG GT TG GG GA AG GC CC CA AA AT TC CC CC CA AA AA AA AA AG GC CC CG GT TC CC CT TC CA AG GT TT TC CG GG GA AT TT TG GC CA AG GT TC CT TG GC CA AA AC CT TC CG GA AC CT TG GC CA AT TG GA AA AG GG GT TG GG GA AA AT TC CG GC CT TA AG GT TA AA AT TC CG GC CG GG GA AT TC CA AG GC CA AT TG GC CC CG GC CG GG GT TG GA AA AT TA AC CG GT TT TC CC CC CG GG GG GC CC CT TT TG GT TA AC CA AC CA AC CC CG GC CC CC CG GT TC CA AC CA AT TC CA AC CG GA AA AA AG GT TT TA AG GT TT TG GC CT TC CT TT TG GA AA AG GT TC CA AC CC CG GA AG GC CT TA AA AT TC CG GC CA AA AG GG GA AG GG GC CA AG
 GG GT TG GC CC CG GA AC CG GG GA AG GT TG GG GC CT TG GG GC CG GA AT TT TG GG GG GG GT TG
In [8]:
print(reference_db[0].kmer_frequencies(k=k))

{'TC': 62, 'TA': 53, 'CC': 99, 'GC': 130, 'AC': 84, 'AG': 105, 'AT': 51, 'CA': 76, 'GA': 110, 'TT': 48, 'AA': 98, 'CG': 122, 'CT': 78, 'GT': 90, 'GG': 167, 'TG': 104}


This information can be convenient to store in a pandas Series object:

In [9]:
pd.Series(reference_db[0].kmer_frequencies(k=k), name=reference_db[0].metadata['id'])

Out[9]:
AA     98
AC     84
AG    105
AT     51
CA     76
CC     99
CG    122
CT     78
GA    110
GC    130
GG    167
GT     90
TA     53
TC     62
TG    104
TT     48
Name: 272225, dtype: int64

To train our taxonomic classifier, we next need to define a few things. First, at what level of taxonomic specificity do we want to classify our sequences? We should expect to achieve higher accuracy at less specific taxonomic levels such as phylum or class, but these are likely to be less informative biologically than more specific levels such as genus or species. Let's start classifying at the phylum level to keep our task simple, since we're working with a small subset of the reference database here. In Greengenes, phylum is the second level of the taxonomy.

Next, how long should our kmers be? We don't have a good idea of this to start with. The longer our kmers, the more likely they are to be specific to certain taxa, which is good because that will help with classification. However, if they get too long it becomes less likely that we'll observe those kmers in sequences that aren't represented in our database because the longer the sequence is the more likely we are to see variation across other organisms that are assigned to the same taxonomy. Based on some of my own work in this area, I'll start us out with 7-mers (i.e., kmers of length 7).

Finally, we'll need to know the value of W, defined above as the set of all possible kmers given our alphabet and the value of k.

As an exercise, I recommend exploring the impact of the value of k and taxonomic_level on the accuracy of our classifier after reading this chapter.

In [10]:
taxonomic_level = 2
k = 7
alphabet = skbio.DNA.nondegenerate_chars


Next, we'll compute a table of the per-sequence kmer counts for all kmers in W for all sequences in our reference database. We'll also store the taxonomic label of each of our reference sequences at our specified taxonomic level. We can store this information in a pandas DataFrame, and then view the first 25 rows of that table.

In [11]:
def get_taxon_at_level(taxon, level):
taxon = [l.strip() for l in taxon.split(';')]
return '; '.join(taxon[:level])

W = compute_W(alphabet, k)

per_sequence_kmer_counts = []
for reference_sequence in reference_db:
kmer_counts = dict.fromkeys(W, 0)
kmer_counts.update(reference_sequence.kmer_frequencies(k=k))
per_sequence_kmer_counts.append(pd.Series(kmer_counts, name=taxon))

per_sequence_kmer_counts = pd.DataFrame(data=per_sequence_kmer_counts).fillna(0).T
per_sequence_kmer_counts[:25]

Out[11]:
k__Bacteria; p__Acidobacteria k__Bacteria; p__Bacteroidetes k__Archaea; p__Crenarchaeota k__Bacteria; p__Bacteroidetes k__Bacteria; p__Bacteroidetes k__Archaea; p__Euryarchaeota k__Bacteria; p__Firmicutes k__Bacteria; p__Proteobacteria k__Bacteria; p__Actinobacteria k__Bacteria; p__Proteobacteria ... k__Bacteria; p__ k__Bacteria; p__Bacteroidetes k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria k__Bacteria; p__Planctomycetes k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria
AAAAAAA 0 0 1 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAAAC 0 0 0 0 0 0 0 0 0 0 ... 0 0 1 0 0 0 0 0 0 0
AAAAAAG 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAAAT 0 0 1 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAACA 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 1 0 0
AAAAACC 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAACG 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAACT 0 0 0 0 0 0 1 0 0 0 ... 0 0 1 0 0 0 0 0 0 0
AAAAAGA 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAAGC 1 1 0 0 0 0 0 0 1 1 ... 1 0 0 0 1 0 0 0 0 0
AAAAAGG 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAAGT 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAATA 0 0 1 0 0 0 1 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAATC 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAATG 0 0 1 0 0 0 1 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAAATT 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAACAA 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 1 0 0
AAAACAC 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAACAG 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 1 0 0 0
AAAACAT 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAACCA 0 0 1 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAACCC 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 1 0 0 0 0
AAAACCG 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAACCT 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
AAAACGA 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0

25 rows × 500 columns

With this information, we'll next compute our "kmer probability table" (EXISTING NAME FOR THIS?). The content of this table will be the probability of observing each kmer in W given a taxon. This is computed based on a few values:

$N$ : The total number of sequences in the training set.

$n(w_i)$ : The number of total sequences containing kmer i.

$P_i$ : The probability of observing kmer i. Initially it might seem as though this would be computed as $n(w_i) / N$, but this neglects the possibility of that a kmer observed in a query sequence might not be represented in our reference database, so a small pseudocount is added to the numerator and denomenator.

$P(w_i | taxon)$ : The probability of observing a kmer given a taxon. Again, it would seem that this would be computed as the proportion of sequences in the taxon containing the kmer, but this would neglect that we'll likely observe kmers in our query sequences that are not represented in our reference database. As pseudocount is therefore added again to the numerator and denominator. This time the pseudocount in the numerator is scaled by how frequent the kmer is in the reference database as a whole: specifically, it is $P_i$.

Our "kmer probability table" is $P(w_i | taxon)$ computed for all kmers in W and all taxa represented in our reference database. We'll compute that and again look at the first 25 rows.

In [12]:
def compute_kmer_probability_table(per_sequence_kmer_counts):
N = len(per_sequence_kmer_counts) # number of training sequences

# number of sequences containing kmer wi
n_wi = per_sequence_kmer_counts.astype(bool).sum(axis=1)
n_wi.name = 'n(w_i)'

# probabilities of observing each kmer
Pi = (n_wi + 0.5) / (N + 1)
Pi.name = 'P_i'

# number of times each taxon appears in training set
taxon_counts = collections.Counter(per_sequence_kmer_counts.columns)
n_taxon_members_containing_kmer = per_sequence_kmer_counts.astype(bool).groupby(level=0, axis=1).sum()

# probabilities of observing each kmer in each taxon
p_wi_t = []
for taxon, count in taxon_counts.items():
p_wi_t.append(pd.Series((n_taxon_members_containing_kmer[taxon] + Pi) / (count + 1), name=taxon))

return pd.DataFrame(p_wi_t).T

In [13]:
kmer_probability_table = compute_kmer_probability_table(per_sequence_kmer_counts)

In [14]:
kmer_probability_table[:25]

Out[14]:
k__Bacteria; p__WS6 k__Bacteria; p__Verrucomicrobia k__Bacteria; p__SAR406 k__Bacteria; p__WS4 k__Bacteria; p__OP8 k__Bacteria; p__Nitrospirae k__Bacteria; p__NKB19 k__Archaea; p__Euryarchaeota k__Bacteria; p__Chlamydiae k__Bacteria; p__Hyd24-12 ... k__Bacteria; p__WS3 k__Bacteria; p__GAL15 k__Bacteria; p__OP11 k__Bacteria; p__AD3 k__Archaea; p__Crenarchaeota k__Bacteria; p__Proteobacteria k__Bacteria; p__Bacteroidetes k__Bacteria; p__Cyanobacteria k__Bacteria; p__ k__Bacteria; p__Spirochaetes
AAAAAAA 0.000254 0.000254 0.000254 0.000381 0.000381 0.000254 0.000381 0.000191 0.000381 0.000381 ... 0.000254 0.000381 0.000381 0.000381 0.125095 0.007098 0.014504 0.153905 0.000381 0.000069
AAAAAAC 0.000783 0.000783 0.667450 0.001175 0.001175 0.334117 0.001175 0.000587 0.001175 0.001175 ... 0.000783 0.001175 0.001175 0.001175 0.000294 0.127676 0.072498 0.000181 0.001175 0.091123
AAAAAAG 0.000519 0.000519 0.000519 0.000778 0.000778 0.000519 0.000778 0.000389 0.000778 0.000778 ... 0.000519 0.000778 0.000778 0.000778 0.000195 0.028380 0.057994 0.077043 0.000778 0.000141
AAAAAAT 0.000417 0.000417 0.000417 0.000626 0.000626 0.000417 0.000626 0.000313 0.000626 0.000626 ... 0.000417 0.000626 0.000626 0.000626 0.125156 0.035470 0.029004 0.230865 0.000626 0.000114
AAAAACA 0.000376 0.000376 0.000376 0.000565 0.000565 0.000376 0.000565 0.000282 0.000565 0.000565 ... 0.000376 0.000565 0.000565 0.000565 0.000141 0.028377 0.043495 0.000087 0.000565 0.000103
AAAAACC 0.001271 0.001271 0.334605 0.001907 0.001907 0.334605 0.001907 0.000954 0.001907 0.001907 ... 0.334605 0.001907 0.001907 0.001907 0.000477 0.226977 0.115997 0.000293 0.001907 0.091256
AAAAACG 0.000295 0.000295 0.000295 0.000442 0.000442 0.333628 0.000442 0.000221 0.000442 0.000442 ... 0.000295 0.000442 0.000442 0.000442 0.000111 0.000006 0.028998 0.000068 0.000442 0.000080
AAAAACT 0.000966 0.000966 0.334300 0.001449 0.001449 0.000966 0.001449 0.000725 0.001449 0.001449 ... 0.000966 0.001449 0.001449 0.001449 0.000362 0.049666 0.029028 0.000223 0.001449 0.000264
AAAAAGA 0.000275 0.000275 0.000275 0.000412 0.000412 0.000275 0.000412 0.000206 0.000412 0.000412 ... 0.000275 0.000412 0.000412 0.000412 0.000103 0.014190 0.043490 0.076986 0.000412 0.000075
AAAAAGC 0.001638 0.001638 0.001638 0.002457 0.002457 0.001638 0.502457 0.001228 0.002457 0.502457 ... 0.334971 0.002457 0.002457 0.002457 0.125614 0.177340 0.130506 0.000378 0.502457 0.091356
AAAAAGG 0.000295 0.000295 0.000295 0.000442 0.000442 0.000295 0.000442 0.000221 0.000442 0.000442 ... 0.000295 0.000442 0.000442 0.000442 0.000111 0.014191 0.072477 0.000068 0.000442 0.000080
AAAAAGT 0.000478 0.000478 0.000478 0.000717 0.000717 0.000478 0.000717 0.000359 0.000717 0.000717 ... 0.000478 0.000717 0.000717 0.000717 0.000179 0.085117 0.014514 0.000110 0.000717 0.000130
AAAAATA 0.000702 0.000702 0.000702 0.001053 0.001053 0.000702 0.001053 0.000526 0.001053 0.001053 ... 0.000702 0.001053 0.001053 0.001053 0.125263 0.000015 0.000031 0.077085 0.001053 0.000191
AAAAATC 0.000193 0.000193 0.333527 0.000290 0.000290 0.000193 0.000290 0.000145 0.000290 0.000290 ... 0.000193 0.000290 0.000290 0.000290 0.000072 0.014189 0.000008 0.000045 0.000290 0.000053
AAAAATG 0.000376 0.000376 0.000376 0.000565 0.000565 0.000376 0.000565 0.000282 0.000565 0.000565 ... 0.000376 0.000565 0.000565 0.000565 0.125141 0.021285 0.014509 0.153933 0.000565 0.000103
AAAAATT 0.000254 0.000254 0.000254 0.000381 0.000381 0.000254 0.000381 0.000191 0.000381 0.000381 ... 0.000254 0.000381 0.000381 0.000381 0.000095 0.028374 0.043489 0.000059 0.000381 0.000069
AAAACAA 0.000254 0.000254 0.000254 0.000381 0.000381 0.000254 0.000381 0.000191 0.000381 0.000381 ... 0.000254 0.000381 0.000381 0.000381 0.000095 0.007098 0.028997 0.000059 0.000381 0.000069
AAAACAC 0.000153 0.000153 0.000153 0.000229 0.000229 0.000153 0.000229 0.000114 0.000229 0.000229 ... 0.000153 0.000229 0.000229 0.000229 0.000057 0.014188 0.014499 0.000035 0.000229 0.000042
AAAACAG 0.000376 0.000376 0.000376 0.000565 0.000565 0.000376 0.000565 0.000282 0.000565 0.000565 ... 0.000376 0.000565 0.000565 0.000565 0.000141 0.035469 0.057987 0.000087 0.000565 0.000103
AAAACAT 0.000153 0.000153 0.000153 0.500229 0.000229 0.000153 0.000229 0.000114 0.000229 0.000229 ... 0.000153 0.000229 0.000229 0.000229 0.000057 0.014188 0.028992 0.000035 0.000229 0.000042
AAAACCA 0.000437 0.000437 0.000437 0.000656 0.000656 0.000437 0.000656 0.000328 0.000656 0.000656 ... 0.000437 0.000656 0.000656 0.000656 0.500164 0.021286 0.043497 0.000101 0.000656 0.000119
AAAACCC 0.000661 0.000661 0.000661 0.000992 0.000992 0.000661 0.000992 0.000496 0.000992 0.000992 ... 0.000661 0.000992 0.000992 0.000992 0.000248 0.007106 0.058000 0.153999 0.000992 0.000180
AAAACCG 0.668040 0.334707 0.334707 0.002060 0.502060 0.334707 0.002060 0.001030 0.002060 0.002060 ... 0.334707 0.002060 0.002060 0.002060 0.000515 0.177334 0.144987 0.000317 0.002060 0.182193
AAAACCT 0.000865 0.000865 0.000865 0.001297 0.001297 0.000865 0.001297 0.000648 0.001297 0.001297 ... 0.000865 0.001297 0.001297 0.001297 0.000324 0.156047 0.159458 0.000200 0.001297 0.000236
AAAACGA 0.000112 0.000112 0.000112 0.000168 0.000168 0.000112 0.000168 0.000084 0.000168 0.000168 ... 0.000112 0.000168 0.000168 0.000168 0.000042 0.000002 0.014498 0.000026 0.000168 0.000031

25 rows × 42 columns

With our kmer probability table we are now ready to classify unknown sequences. We'll begin by defining some query sequences. We'll pull these at random from our reference sequences, which means that some of the query sequences will be represented in our reference database and some won't be. This is the sitatuation that is typically encountered in practice. To simulate real-world 16S taxonomy classification tasks, we'll also trim out 200 bases of our reference sequences since (as of this writing) we typically don't obtain full-length 16S sequences from a DNA sequencing instrument.

In [15]:
queries = []
for e in skbio.io.read(qdr.get_reference_sequences(), format='fasta', constructor=skbio.DNA):
if e.has_degenerates():
# For the purpose of this lesson, we're going to ignore sequences that contain
# degenerate characters (i.e., characters other than A, C, G, or T)
continue
e = e[100:300]
queries.append(e)
# can't figure out why np.random.choice isn't working here...
np.random.shuffle(queries)
queries = queries[:50]

In [16]:
queries[0]

Out[16]:
DNA
---------------------------------------------------------------------
'description': ''
'id': '36764'
Stats:
length: 200
has gaps: False
has degenerates: False
has non-degenerates: True
GC-content: 55.50%
---------------------------------------------------------------------
0   AGAATCTGCC TTCAGGATGG GGACAACAGT GGGAAACTGC TGCTAATCCC CAATAAGCCG
60  AGAGGTGAAA GAGTTTCGCC TGAAGAGGAG CTCGCGTCCG ATTAGCTAGT TGGTGGGGTA
120 AAGGCTTACC AAGGCGACGA TCGGTAGCTG GTCTGAGAGG ACGATCAGCC ACACTGGGAC
180 TGAGACACGG CCCAGACTCC

For a given query sequence, its taxonomy will be classified as follows. First, the set of all kmers will be extracted from the sequence. This is referred to as $V$. Then, for all taxa in the kmer probability table, the probability of observing the query sequence will be computed given that taxon: $P(query | taxon)$. This is computed as the product of all its kmer probabilities for the given taxon. (It should be clear based on this formula why it was necessary to add pseudocounts when computing our kmer probability table - if not, kmer probabilities of zero would result in a zero probability of the sequence being derived from that taxon at this step.)

After computing $P(query | taxon)$ for all taxa, the taxonomy assignment return is simply the one achieving the maximum probability. Here we'll classify a sequence and look at the resulting taxonomy assignment.

In [17]:
def classify_V(V, kmer_probability_table):
P_S_t = [] # probability of the sequence given the taxon
for taxon in kmer_probability_table:
kmer_probabilities = kmer_probability_table[taxon]
probability = 1.0
for v_i in V:
probability *= kmer_probabilities[v_i]
P_S_t.append((probability, taxon))
return max(P_S_t)[1], V

def classify_sequence(query_sequence, kmer_probability_table, k):
V = list(map(str, query_sequence.iter_kmers(k=k)))
return classify_V(V, kmer_probability_table)

In [18]:
taxon_assignment, V = classify_sequence(queries[0], kmer_probability_table, k)
print(taxon_assignment)

k__Bacteria; p__Proteobacteria


Since we know the actual taxonomy assignment for this sequence, we can look that up in our reference database. Was your assignment correct? Try this with a few query sequences and keep track of how many times the classifier achieved the correct assignment.

In [19]:
get_taxon_at_level(reference_taxonomy[queries[0].metadata['id']], taxonomic_level)

Out[19]:
'k__Bacteria; p__Cyanobacteria'

Because the query and reference sequences that were working with were randomly selected from the full reference database, each time you run this notebook you should observe different results. Chances are however that if you run the above steps multiple times you'll get the wrong taxonomy assignment at least some of the time. Up to this point, we've left out an important piece of information: how confident should we be in our assignment, or in other workds, how dependent is our taxonomy assignment on our specific query? If there were slight differences in our query (e.g., because we observed a very closely related organism, such as one of the same species but a different strain, or because we sequenced a different region of the 16S sequence) would we obtain the same taxonomy assignment? If so, we should have higher confidence in our assignment. If not, we should have lower confidence in our assignment.

We can quantify confidence using an approach called bootstrapping. With a bootstrap approach, we'll get our taxonomy assignment as we did above, but then for some user-specified number of times, we'll create random subsets of V sampled with replacement (DEFINE THIS). We'll then assign taxonomy each random subset of V, and count the number of times the resulting taxonomy assignment is the same that we achieved when assigning taxonomy to V. The count divided by the number of iterations we've chosen to run will be our confidence value. If the assignments are often the same we'll have a high confidence value. If the assignments are often different, we'll have a low confidence value.

Let's now assign taxonomy and compute a confidence for that assignment.

In [20]:
def classify_sequence_with_confidence(sequence, kmer_probability_table, k,
confidence_iterations=100):
taxon, V = classify_sequence(sequence, kmer_probability_table, k)

count_same_taxon = 0
subsample_size = int(len(V) * 0.1)
for i in range(confidence_iterations):
subsample_V = np.random.choice(V, subsample_size, replace=True)
subsample_taxon, _ = classify_V(subsample_V, kmer_probability_table)
if taxon == subsample_taxon:
count_same_taxon += 1
confidence = count_same_taxon / confidence_iterations

return (taxon, confidence)

In [21]:
taxon_assignment, confidence = classify_sequence_with_confidence(queries[0], kmer_probability_table, k)
print(taxon_assignment)
print(confidence)

k__Bacteria; p__Proteobacteria
0.44


How did the computed confidence compare to the accuracy taxonomy assignment?

We don't have an a priori idea for what good versus bad confidence scores are, but we can use our reference database to explore that. We might want this information so we can come up with a confidence threshold, above which we would accept a taxonomy assignment and below which we might reject it. To explore this, let's compute taxonomy assignments and confidence for all of our query sequences and then see what the distributions of confidence scores look like for correct assignments and incorrect assignments.

In [22]:
correct_assignment_confidences = []
incorrect_assignment_confidences = []
summary = []

for query in queries:
predicted_taxonomy, confidence = classify_sequence_with_confidence(query, kmer_probability_table, k)
if actual_taxonomy == predicted_taxonomy:
correct_assignment_confidences.append(confidence)
else:
incorrect_assignment_confidences.append(confidence)

summary.append([predicted_taxonomy, actual_taxonomy, confidence])
summary = pd.DataFrame(summary, columns=['Predicted taxonomy', 'Actual taxonomy', 'Confidence'])

In [23]:
import seaborn as sns

ax = sns.boxplot(data=[correct_assignment_confidences, incorrect_assignment_confidences])
ax = sns.swarmplot(data=[correct_assignment_confidences, incorrect_assignment_confidences], color="black")
_ = ax.set_xticklabels(['Correct assignments', 'Incorrect assignments'])
_ = ax.set_ylabel('Confidence')


What does this plot tell you about how well setting a confidence threshold is likely to work? If you never wanted to reject a correct assignment, how often would you accept an incorrect assignment? If you never wanted to accept an incorrect assignment, how often would you reject a correct assignment?

In [24]:
summary # maybe explore whether certain taxa are more frequently wrong than others...

Out[24]:
Predicted taxonomy Actual taxonomy Confidence
0 k__Bacteria; p__Proteobacteria k__Bacteria; p__Cyanobacteria 0.45
1 k__Bacteria; p__Proteobacteria k__Archaea; p__Euryarchaeota 0.34
2 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.89
3 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.92
4 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.92
5 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.95
6 k__Bacteria; p__Bacteroidetes k__Bacteria; p__Bacteroidetes 0.66
7 k__Bacteria; p__Planctomycetes k__Bacteria; p__Planctomycetes 0.51
8 k__Bacteria; p__Planctomycetes k__Bacteria; p__Planctomycetes 0.92
9 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.78
10 k__Bacteria; p__Proteobacteria k__Bacteria; p__Acidobacteria 0.30
11 k__Bacteria; p__Firmicutes k__Bacteria; p__Acidobacteria 0.48
12 k__Bacteria; p__Bacteroidetes k__Bacteria; p__Bacteroidetes 0.80
13 k__Bacteria; p__Proteobacteria k__Archaea; p__Euryarchaeota 0.67
14 k__Bacteria; p__Proteobacteria k__Bacteria; p__SAR406 0.60
15 k__Bacteria; p__Proteobacteria k__Bacteria; p__Verrucomicrobia 0.51
16 k__Bacteria; p__Proteobacteria k__Bacteria; p__Verrucomicrobia 0.45
17 k__Bacteria; p__Actinobacteria k__Bacteria; p__Actinobacteria 0.67
18 k__Bacteria; p__Bacteroidetes k__Bacteria; p__Bacteroidetes 0.78
19 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.77
20 k__Bacteria; p__Actinobacteria k__Bacteria; p__OP8 0.57
21 k__Bacteria; p__Bacteroidetes k__Bacteria; p__Bacteroidetes 0.91
22 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.68
23 k__Bacteria; p__Bacteroidetes k__Bacteria; p__Bacteroidetes 0.94
24 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 1.00
25 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.87
26 k__Bacteria; p__Firmicutes k__Bacteria; p__NC10 0.49
27 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.91
28 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 1.00
29 k__Bacteria; p__Proteobacteria k__Bacteria; p__Chloroflexi 0.52
30 k__Bacteria; p__Actinobacteria k__Bacteria; p__Actinobacteria 0.44
31 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.98
32 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.90
33 k__Bacteria; p__Proteobacteria k__Bacteria; p__GN04 0.70
34 k__Bacteria; p__Proteobacteria k__Bacteria; p__SAR406 0.57
35 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.99
36 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.92
37 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.87
38 k__Bacteria; p__Bacteroidetes k__Bacteria; p__Bacteroidetes 0.99
39 k__Bacteria; p__Proteobacteria k__Bacteria; p__Chloroflexi 0.55
40 k__Bacteria; p__Acidobacteria k__Bacteria; p__Acidobacteria 0.36
41 k__Bacteria; p__Proteobacteria k__Bacteria; p__Acidobacteria 0.37
42 k__Bacteria; p__Planctomycetes k__Bacteria; p__Planctomycetes 0.47
43 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.92
44 k__Bacteria; p__Actinobacteria k__Bacteria; p__Actinobacteria 0.56
45 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 1.00
46 k__Bacteria; p__Firmicutes k__Bacteria; p__Firmicutes 0.92
47 k__Bacteria; p__Proteobacteria k__Bacteria; p__Chloroflexi 0.63
48 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 1.00
49 k__Bacteria; p__Proteobacteria k__Bacteria; p__Proteobacteria 0.52