In [1]:
%matplotlib inline
from os.path import join, exists, expandvars
import pandas as pd
from IPython.display import display
from tax_credit.plotting_functions import (pointplot_from_data_frame,
boxplot_from_data_frame,
heatmap_from_data_frame,
per_level_kruskal_wallis,
beta_diversity_pcoa,
average_distance_boxplots,
rank_optimized_method_performance_by_dataset)
from tax_credit.eval_framework import (evaluate_results,
method_by_dataset_a1,
method_by_dataset_a2,
parameter_comparisons,
merge_expected_and_observed_tables)
This is the only cell that you will need to edit to generate basic reports locally. After editing this cell, you can run all cells in this notebook to generate your analysis report. This will take a few minutes to run, as results are computed at multiple taxonomic levels.
Values in this cell will not need to be changed, with the exception of project_dir, to generate the default results contained within tax-credit. To analyze results separately from the tax-credit precomputed results, other variables in this cell will need to be set.
In [2]:
## project_dir should be the directory where you've downloaded (or cloned) the
## tax-credit repository.
project_dir = expandvars("../../")
## expected_results_dir contains expected composition data in the structure
## expected_results_dir/<dataset name>/<reference name>/expected/
expected_results_dir = join(project_dir, "data/precomputed-results/simulated-community")
## simulated_results_fp designates the files to which summary results are written.
## If this file exists, it can be read in to generate results plots, instead
## of computing new scores.
simulated_results_fp = join(expected_results_dir, 'simulated_results.tsv')
## results_dirs should contain the directory or directories where
## results can be found. By default, this is the same location as expected.
results_dirs = [expected_results_dir]
## Minimum number of times an OTU must be observed for it to be included in analyses. Edit this
## to analyze the effect of the minimum count on taxonomic results.
min_count = 1
## Define the range of taxonomic levels over which to compute accuracy scores.
## The default given below will compute order (level 2) through species (level 6)
taxonomy_level_range = range(2,7)
Next we'll use the paths defined above to find all of the tables that will be compared. These include the pre-computed result tables (i.e., the ones that the new methods will be compared to), the expected result tables (i.e., the tables containing the known composition of the simulated microbial communities), and the query result tables (i.e., the tables generated with the new method(s) that we want to compare to the pre-computed result tables).
Note: if you have added additional methods and are attempting to recompute results, set force=True.
This cell will take a few minutes to run, so hold onto your hat.
In [3]:
mock_results = evaluate_results(results_dirs,
expected_results_dir,
simulated_results_fp,
taxonomy_level_range=range(2,7),
min_count=min_count,
taxa_to_keep=None,
md_key='taxonomy',
subsample=False,
filename_pattern='expected-composition.biom',
force=False)
Restrict analyses to a set of datasets or references: e.g., exclude taxonomy assignments made for purpose of reference database comparisons.
In [4]:
mock_results = mock_results[mock_results['Reference'].isin(['gg_13_8_otus', 'unite_20.11.2016'])]
In this evaluation, we compute and summarize precision, recall, and F-measure of each result (pre-computed and query) based on the known composition of the mock communities. We then summarize the results in two ways: first with boxplots, and second with a table of the top methods based on their F-measures.
This is a qualitative evaluation, effectively telling us about the ability of the different methods to report the taxa that are present in each sample. These metrics are not concerned with the abundance of the different taxa.
As a first step, we will evaluate how accuracy metrics perform at each taxonomic level for each method within each reference dataset type.
In [5]:
color_pallette={'expected': 'red', 'rdp': 'seagreen', 'sortmerna': 'gray', 'uclust': 'blue', 'blast': 'purple'}
y_vars = ["Precision", "Recall", "F-measure"]
In [87]:
pointplot_from_data_frame(mock_results, "Level", y_vars,
group_by="Reference", color_by="Method",
color_pallette=color_pallette)
In [88]:
result = per_level_kruskal_wallis(mock_results, y_vars, group_by='Method',
dataset_col='Reference', level_name='Level',
levelrange=range(2,7), alpha=0.05,
pval_correction='fdr_bh')
result
Out[88]:
Principal coordinate analysis offers a neat way to assess the relative performance of multiple methods to reconstruct expected compositions. Methods that cluster with the "expected" composition probably outperform those that appear more distant on a PCoA plot. First, we need to merge biom tables from each method/parameter configuration for each dataset/reference/level combination, so that we can compare each method/parameter as a separate "sample".
Note: if you have added additional methods and are attempting to recompute results, set force=True.
In [89]:
merge_expected_and_observed_tables(expected_results_dir, results_dirs, taxonomy_level=6,
filename_pattern='expected-composition.biom', force=False)
Now we can manually select which table we want to view. This will output a Bray-Curtis PCoA plot, in addition to ANOSIM test results, which indicate whether at least two methods are significantly different from each other.
These plots are useful for visualizing the relative performance of different methods and their configurations relative to each other and to expected compositions, but are primarily a qualitative technique and do not really give us an idea of whether method X actually performs better than method Y.
Note that 2D plots will only appear if you are running notebooks locally. If viewing static notebooks online, make sure you are viewing this notebook in nbviewer. (if viewing on GitHub, just copy the URL and paste into the search bar in nbviewer.)
In [6]:
table = join(expected_results_dir, 'sake', 'gg_13_8_otus', 'merged_table.biom')
sample_md, results, pc, dm = beta_diversity_pcoa(table, method="braycurtis", dim=2,
permutations=99, col='method',
colormap=color_pallette)
In [7]:
table = join(expected_results_dir, 'wine', 'unite_20.11.2016', 'merged_table.biom')
sample_md, results, pc, dm = beta_diversity_pcoa(table, method="braycurtis", dim=2,
permutations=99, col='method',
colormap=color_pallette)
Let's generate more Bray-Curtis distance plots at genus level. Unlike the species-level plots, in which samples cluster primarily by method configuration, and some method configurations perform considerably better than others, we see at genus level that samples cluster instead according to sample type with less regard to method configuration, indicating that most method configurations perform well at genus level, but few perform well at species level.
In [91]:
merge_expected_and_observed_tables(expected_results_dir, results_dirs, taxonomy_level=5,
filename_pattern='expected-composition.biom', force=False,
biom_fp='merged_table_L5.biom')
In [8]:
table = join(expected_results_dir, 'sake', 'gg_13_8_otus', 'merged_table_L5.biom')
sample_md, results, pc, dm = beta_diversity_pcoa(table, method="braycurtis", dim=2,
permutations=99, col='method',
colormap=color_pallette)
In [9]:
table = join(expected_results_dir, 'wine', 'unite_20.11.2016', 'merged_table_L5.biom')
sample_md, results, pc, dm = beta_diversity_pcoa(table, method="braycurtis", dim=2,
permutations=99, col='method',
colormap=color_pallette)
Generates violin plots of distribution of distances between expected composition and predicted compositions for each method (all parameter configurations) across all samples/datasets, and pairwise Mann Whitney U tests between these distributions.
In [64]:
results = average_distance_boxplots(expected_results_dir, paired=False, use_best=False)
results
Out[64]:
Reports the top-performing parameter configuration for each method, violin plots of distribution of distances between expected composition and predicted compositions for the top parameter for each method across all samples/datasets, and pairwise paired Wilcoxon signed rank tests between these distributions.
In [65]:
results = average_distance_boxplots(expected_results_dir)
results
Out[65]:
In [68]:
mock_results_6 = mock_results[mock_results['Level'] == 6]
In [73]:
for i in mock_results_6['Dataset'].unique():
print(i)
best = method_by_dataset_a1(mock_results_6, i)
display(best)
Now we can determine which parameter configuration performed best for each method. Count best values in each column indicate how many samples a given method achieved within one mean absolute deviation of the best result (which is why they may sum to more than the total number of samples).
In [74]:
for method in mock_results_6['Method'].unique():
top_params = parameter_comparisons(mock_results_6, method)
print(method)
display(top_params[:5])