{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# [2.6]() Machine learning in bioinformatics (work-in-progress) [edit]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Table of Contents**\n",
"0. [Defining a classification problem](#1)\n",
" 0. [Naive Bayes classifiers](#1.1)\n",
" 0. [Random Forest classifiers](#1.2)\n",
"0. [Defining a dimensionality reduction problem](#2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**This chapter is currently a work-in-progress, and is incomplete.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Machine learning algorithms are commonly used in bioinformatics for a variety of tasks. Typically, the common thread in these tasks is that the user would like the algorithm to assist in the identification of patterns in a complex data set. In this chapter we'll implement a few machine learning algorithms 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](http://scikit-learn.org/) is a popular and well-documented Python library for machine learning which many bioinformatics researchers and software developers use in their work."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These algorithms generally work beginning with a collection of samples and some user-defined features of those samples. These data are typically represented in a matrix, where samples are the rows and features are the columns. There are a few different high-level tasks that are common in machine learning, including classification, regression, and dimensionality reduction. In a classification task, a user provides examples of data that fall into certain discrete classes (for example, _healthy_ and _disease_), and tries to have the computer develop a model that can differentiate those classes based on the defined features. If successful, the resulting model could be applied to data where the class isn't known ahead of time, in attempt to predict the class from the features. A regression task is similar, except that a continuous value will be predicted rather than a discrete value. Dimensionality reduction tasks, on the other hand, generally don't have classes or labels assigned ahead of time, and the user is hoping to identify which samples are most similar to each other based on new features that are defined by the algorithm. The goal here might be to reduce the number of features from thousands or more to around two or three that explain most of the variation in the data. This allows the user to explore the samples visually, for example in a scatter plot, which would not be feasible if there were thousands of features."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this chapter we'll explore two classification algorithms and one dimensionality reduction task in the context of some real-world examples."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## [2.6.1](#1) Defining a classification problem [edit]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Table of Contents**\n",
"0. [Naive Bayes classifiers](#1.1)\n",
"0. [Random Forest classifiers](#1.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"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](2.ipynb), so it is likely worth spending a few minutes skimming that chapter if it's not fresh in your mind."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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](http://greengenes.secondgenome.com/) database, which we'll access using [QIIME default reference project](https://github.com/biocore/qiime-default-reference). 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](2.ipynb).)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This time, instead of using sequence alignment to identify the most likely taxonomic origin of a sequence, we'll train classifiers by building [kmer](../5/2.ipynb#2)-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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's jump in..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### [2.6.1.1](#1.1) Naive Bayes classifiers [edit]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"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](../5/2.ipynb#2)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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](2.ipynb#2)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%pylab inline\n",
"\n",
"from IPython.core import page\n",
"page.page = print\n",
"\n",
"import pandas as pd\n",
"import skbio\n",
"import numpy as np\n",
"import itertools\n",
"import collections"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from iab.algorithms import load_taxonomy_reference_database\n",
"\n",
"%psource load_taxonomy_reference_database"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"reference_taxonomy, reference_db = load_taxonomy_reference_database()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"reference_db[0]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"reference_db[-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll select a random subset of the reference database to work with here."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"reference_db = np.random.choice(reference_db, 500, replace=False)\n",
"print(\"%s sequences are present in the subsampled database.\" % len(reference_db))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)?"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"alphabet = skbio.DNA.nondegenerate_chars\n",
"k = 2\n",
"\n",
"def compute_W(alphabet, k):\n",
" return set(map(''.join, itertools.product(alphabet, repeat=k)))\n",
"\n",
"W = compute_W(alphabet, k)\n",
"print('Alphabet contains the characters: %s' % ', '.join(alphabet))\n",
"print('For an alphabet size of %d, W contains %d length-%d kmers.' % (len(alphabet), len(W), k))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"kmers = reference_db[0].iter_kmers(k=k)\n",
"for kmer in kmers:\n",
" print(kmer, end=' ')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"print(reference_db[0].kmer_frequencies(k=k))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This information can be convenient to store in a pandas ``Series`` object:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"pd.Series(reference_db[0].kmer_frequencies(k=k), name=reference_db[0].metadata['id'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"taxonomic_level = 2\n",
"k = 7\n",
"alphabet = skbio.DNA.nondegenerate_chars"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def get_taxon_at_level(taxon, level):\n",
" taxon = [l.strip() for l in taxon.split(';')]\n",
" return '; '.join(taxon[:level])\n",
"\n",
"W = compute_W(alphabet, k)\n",
"\n",
"per_sequence_kmer_counts = []\n",
"for reference_sequence in reference_db:\n",
" taxon = get_taxon_at_level(reference_sequence.metadata['taxonomy'], taxonomic_level)\n",
" kmer_counts = dict.fromkeys(W, 0)\n",
" kmer_counts.update(reference_sequence.kmer_frequencies(k=k))\n",
" per_sequence_kmer_counts.append(pd.Series(kmer_counts, name=taxon))\n",
"\n",
"per_sequence_kmer_counts = pd.DataFrame(data=per_sequence_kmer_counts).fillna(0).T\n",
"per_sequence_kmer_counts[:25]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$N$ : The total number of sequences in the training set."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$n(w_i)$ : The number of total sequences containing kmer _i_."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$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$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def compute_kmer_probability_table(per_sequence_kmer_counts):\n",
" N = len(per_sequence_kmer_counts) # number of training sequences\n",
"\n",
" # number of sequences containing kmer wi\n",
" n_wi = per_sequence_kmer_counts.astype(bool).sum(axis=1)\n",
" n_wi.name = 'n(w_i)'\n",
"\n",
" # probabilities of observing each kmer\n",
" Pi = (n_wi + 0.5) / (N + 1)\n",
" Pi.name = 'P_i'\n",
"\n",
" # number of times each taxon appears in training set\n",
" taxon_counts = collections.Counter(per_sequence_kmer_counts.columns)\n",
" n_taxon_members_containing_kmer = per_sequence_kmer_counts.astype(bool).groupby(level=0, axis=1).sum()\n",
"\n",
" # probabilities of observing each kmer in each taxon\n",
" p_wi_t = []\n",
" for taxon, count in taxon_counts.items():\n",
" p_wi_t.append(pd.Series((n_taxon_members_containing_kmer[taxon] + Pi) / (count + 1), name=taxon))\n",
"\n",
" return pd.DataFrame(p_wi_t).T"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"kmer_probability_table = compute_kmer_probability_table(per_sequence_kmer_counts)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"kmer_probability_table[:25]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"from iab.algorithms import load_taxonomy_query_sequences\n",
"\n",
"%psource load_taxonomy_query_sequences"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"import random\n",
"\n",
"queries = load_taxonomy_query_sequences()\n",
"queries = random.sample(queries, k=50)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"queries[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def classify_V(V, kmer_probability_table):\n",
" P_S_t = [] # probability of the sequence given the taxon\n",
" for taxon in kmer_probability_table:\n",
" kmer_probabilities = kmer_probability_table[taxon]\n",
" probability = 1.0\n",
" for v_i in V:\n",
" probability *= kmer_probabilities[v_i]\n",
" P_S_t.append((probability, taxon))\n",
" return max(P_S_t)[1], V\n",
"\n",
"def classify_sequence(query_sequence, kmer_probability_table, k):\n",
" V = list(map(str, query_sequence.iter_kmers(k=k)))\n",
" return classify_V(V, kmer_probability_table)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"taxon_assignment, V = classify_sequence(queries[0], kmer_probability_table, k)\n",
"print(taxon_assignment)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"get_taxon_at_level(reference_taxonomy[queries[0].metadata['id']], taxonomic_level)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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 words, 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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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 to 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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's now assign taxonomy and compute a confidence for that assignment."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"def classify_sequence_with_confidence(sequence, kmer_probability_table, k,\n",
" confidence_iterations=100):\n",
" taxon, V = classify_sequence(sequence, kmer_probability_table, k)\n",
"\n",
" count_same_taxon = 0\n",
" subsample_size = int(len(V) * 0.1)\n",
" for i in range(confidence_iterations):\n",
" subsample_V = np.random.choice(V, subsample_size, replace=True)\n",
" subsample_taxon, _ = classify_V(subsample_V, kmer_probability_table)\n",
" if taxon == subsample_taxon:\n",
" count_same_taxon += 1\n",
" confidence = count_same_taxon / confidence_iterations\n",
"\n",
" return (taxon, confidence)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"taxon_assignment, confidence = classify_sequence_with_confidence(queries[0], kmer_probability_table, k)\n",
"print(taxon_assignment)\n",
"print(confidence)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How did the computed confidence compare to the accuracy taxonomy assignment?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"correct_assignment_confidences = []\n",
"incorrect_assignment_confidences = []\n",
"summary = []\n",
"\n",
"for query in queries:\n",
" predicted_taxonomy, confidence = classify_sequence_with_confidence(query, kmer_probability_table, k)\n",
" actual_taxonomy = get_taxon_at_level(reference_taxonomy[query.metadata['id']], taxonomic_level)\n",
" if actual_taxonomy == predicted_taxonomy:\n",
" correct_assignment_confidences.append(confidence)\n",
" else:\n",
" incorrect_assignment_confidences.append(confidence)\n",
"\n",
" summary.append([predicted_taxonomy, actual_taxonomy, confidence])\n",
"summary = pd.DataFrame(summary, columns=['Predicted taxonomy', 'Actual taxonomy', 'Confidence'])"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"import seaborn as sns\n",
"\n",
"ax = sns.boxplot(data=[correct_assignment_confidences, incorrect_assignment_confidences])\n",
"ax = sns.swarmplot(data=[correct_assignment_confidences, incorrect_assignment_confidences], color=\"black\")\n",
"_ = ax.set_xticklabels(['Correct assignments', 'Incorrect assignments'])\n",
"_ = ax.set_ylabel('Confidence')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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?"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"summary # maybe explore whether certain taxa are more frequently wrong than others..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### [2.6.1.2](#1.2) Random Forest classifiers [edit]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Coming soon..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## [2.6.2](#2) Defining a dimensionality reduction problem [edit]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"[This content](../3/1.ipynb#4.3) will be adapted and ported here."
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 4
}