<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Evaluate mock community classification accuracy for nb-extra
The purpose of this notebook is to evaluate taxonomic classification accuracy of mock communities using different classification methods.
Contains an additional section that uses CART analysis.
Prepare the environment
First we'll import various functions that we'll need for generating the report.
Step1:  Configure local environment-specific values
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.
Step2:  Find mock community pre-computed tables, expected tables, and "query" tables
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 mock 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
Step3:  Restrict analyses to a set of datasets or references
Step4:  Compute and summarize precision, recall, and F-measure for mock communities
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
Step5:  CART Analysis
In this section we will use Classification and Regression Trees to try to pick good parameters for the naïve Bayes classifier. In each case, we pick the path to the classification leaf that yields the highest expected F-measure. Also, we unconventionally turn of pruning, so that all parameters are sepecified. This has the effect of picking arbitrary parameters towards the leaves of the decision tree where it doesn't matter as much which parameters we choose.
This section requires the additional dependencies of rpy2 in Python and rpart in R. If you do not wish to install those dependencies, skip the CART Analysis section.
Step6:  Split the Parameter String and Aggregate by Community
Step7:  Kruskal-Wallis between-method accuracy comparisons
Kruskal-Wallis FDR-corrected p-values comparing classification methods at each level of taxonomic assignment
Step8:  Violin plots of per-level accuracy
Heatmaps show the performance of individual method/parameter combinations at each taxonomic level, in each reference database (i.e., for bacterial and fungal mock communities individually).
Now we will focus on results at species level (for genus level, change to level 5)
Step9:  Method Optimization
Which method/parameter configuration performed "best" for a given score? We can rank the top-performing configuration by dataset, method, and taxonomic level.
First, the top-performing method/configuration combination by dataset.
Step10:  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).
Step11:  Optimized method performance
And, finally, which method performed best at each individual taxonomic level for each reference dataset (i.e., for across all fungal and bacterial mock communities combined)?
For this analysis, we rank the top-performing method/parameter combination for each method at family through species levels. Methods are ranked by top F-measure, and the average value for each metric is shown (rather than count best as above). F-measure distributions are plotted for each method, and compared using paired t-tests with FDR-corrected P-values. This cell does not need to be altered, unless if you wish to change the metric used for sorting best methods and for plotting.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from os.path import join, exists, expandvars
import pandas as pd
from IPython.display import display, Markdown
import seaborn.xkcd_rgb as colors
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,
                                       parameter_comparisons,
                                       merge_expected_and_observed_tables,
                                       filter_df)
Explanation: Evaluate mock community classification accuracy for nb-extra
The purpose of this notebook is to evaluate taxonomic classification accuracy of mock communities using different classification methods.
Contains an additional section that uses CART analysis.
Prepare the environment
First we'll import various functions that we'll need for generating the report.
End of explanation
## project_dir should be the directory where you've downloaded (or cloned) the 
## tax-credit repository. 
project_dir = join('..', '..')
## 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/", "mock-community")
## mock_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.
mock_results_fp = join(expected_results_dir, 'broad_sweep_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 included with the project. If other results should be included, 
## absolute paths to those directories should be added to this list.
results_dirs = [expected_results_dir]
## directory containing mock community data, e.g., feature table without taxonomy
mock_dir = join(project_dir, "data", "mock-community")
## 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)
dataset_ids = ['mock-' + str(m) for m in (3, 12, 18, 22, 24, '26-ITS1', '26-ITS9')]
Explanation: Configure local environment-specific values
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.
End of explanation
mock_results = evaluate_results(results_dirs, 
                                expected_results_dir, 
                                mock_results_fp, 
                                mock_dir,
                                taxonomy_level_range=range(2,7), 
                                min_count=min_count,
                                taxa_to_keep=None, 
                                md_key='taxonomy', 
                                subsample=False,
                                per_seq_precision=True,
                                exclude=['other'],
                                method_ids=['nb-extra'],
                                append=False,
                                force=False)
mock_results['Reference'].unique()
Explanation: Find mock community pre-computed tables, expected tables, and "query" tables
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 mock 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 to add, set append=True. If you are attempting to recompute pre-computed results, set force=True.
This cell will take a few minutes to run if new results are being added, so hold onto your hat. If you are attempting to re-compute everything, it may take an hour or so, so go take a nap.
End of explanation
#mock_results = filter_df(mock_results, column_name='Reference',
#                         values=['gg_13_8_otus_amplicon', 'gg_13_8_otus_read', 'gg_13_8_otus_full'],
#                         exclude=False)
mock_results = mock_results.reset_index(drop=True)
Explanation: Restrict analyses to a set of datasets or references: e.g., exclude taxonomy assignments made for purpose of reference database comparisons. This can be performed as shown below — alternatively, specific reference databases, datasets, methods, or parameters can be chosen by setting dataset_ids, reference_ids, method_ids, and parameter_ids in the evaluate_results command above.
End of explanation
color_pallette={
    'nb-extra': 'black'
}
y_vars = ["Precision", "Recall", "F-measure", "Taxon Accuracy Rate", "Taxon Detection Rate"]
pointplot_from_data_frame?
pointplot_from_data_frame(mock_results, "Level", y_vars,
                          "Reference", "Method", color_pallette)
Explanation: Compute and summarize precision, recall, and F-measure for mock communities
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. Higher scores = better accuracy
As a first step, we will evaluate average method performance at each taxonomic level for each method within each reference dataset type.
Note that, as parameter configurations can cause results to vary widely, average results are not a good representation of the "best" results. See here for results using optimized parameters for each method.
First we will define our color palette and the variables we want to plot. Via seaborn, we can apply the xkcd crowdsourced color names. If that still doesn't match your hue, use hex codes.
End of explanation
mock_results['Reference'].unique()
from itertools import product
from pandas import DataFrame, concat, to_numeric
from numpy import mean
import rpy2
%load_ext rpy2.ipython
%R require(rpart)
Explanation: CART Analysis
In this section we will use Classification and Regression Trees to try to pick good parameters for the naïve Bayes classifier. In each case, we pick the path to the classification leaf that yields the highest expected F-measure. Also, we unconventionally turn of pruning, so that all parameters are sepecified. This has the effect of picking arbitrary parameters towards the leaves of the decision tree where it doesn't matter as much which parameters we choose.
This section requires the additional dependencies of rpy2 in Python and rpart in R. If you do not wish to install those dependencies, skip the CART Analysis section.
End of explanation
columns = ['Alpha', 'Class-Prior', 'N-Features', 'Ngram-Range', 'Norm', 'Use-IDF', 'Confidence']
params = DataFrame((s.split(':') for s in mock_results['Parameters']), columns=columns)
keepers = ['Dataset', 'Level', 'Reference']
raw_param_results = concat([mock_results[keepers + ['F-measure']], params], axis=1)
raw_param_results = raw_param_results.apply(to_numeric, errors='ignore')
param_results = raw_param_results.groupby(keepers + columns, as_index=False).mean()
len(param_results)
%%R
recommend_params <- function(data, prior, levels, references)
{
    data = data[data[,"Reference"] %in% references,]
    data = data[data[,"Class.Prior"] == prior,]
    data = data[data[,"Level"] %in% levels,]
    fit <- rpart(F.measure ~ Confidence + Use.IDF + Ngram.Range + N.Features + Alpha + Reference + Norm, 
                 data=data,
                 method="anova",
                 control=rpart.control(cp=0))
    rightmost_leaf <- fit$frame[fit$frame[,"yval"] == max(fit$frame[,"yval"]),]
    path.rpart(fit, as.numeric(rownames(rightmost_leaf)))
}
priors = ['uniform', 'prior']
reference_sets = [
    ['gg_13_8_otus_amplicon', 'gg_13_8_otus_full', 'gg_13_8_otus_read'],
    ['unite_20.11.2016_clean_amplicon', 'unite_20.11.2016_clean_full',
    'unite_20.11.2016_clean_read']
]
level_sets = [[2,3,4,5], [6]]
for prior, levels, references in product(priors, level_sets, reference_sets):
    display(Markdown("Prior: `" + prior + '`'))
    display(Markdown("References: `"  + str(references) + '`'))
    display(Markdown("Levels: `" + str(levels) + '`'))
    %R -i param_results,prior,levels,references recommend_params(param_results, prior, levels, references)
Explanation: Split the Parameter String and Aggregate by Community
End of explanation
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
Explanation: Kruskal-Wallis between-method accuracy comparisons
Kruskal-Wallis FDR-corrected p-values comparing classification methods at each level of taxonomic assignment
End of explanation
mock_results_6 = mock_results[mock_results['Level'] == 6]
boxplot_from_data_frame?
boxplot_from_data_frame(mock_results_6, group_by="Method", metric="Precision", color_palette=color_pallette)
boxplot_from_data_frame(mock_results_6, group_by="Method", metric="Recall", color_palette=color_pallette)
boxplot_from_data_frame(mock_results_6, group_by="Method", metric="F-measure", color_palette=color_pallette)
boxplot_from_data_frame(mock_results_6, group_by="Method", metric="Taxon Accuracy Rate", color_palette=color_pallette)
boxplot_from_data_frame(mock_results_6, group_by="Method", metric="Taxon Detection Rate", color_palette=color_pallette)
Explanation: Violin plots of per-level accuracy
Heatmaps show the performance of individual method/parameter combinations at each taxonomic level, in each reference database (i.e., for bacterial and fungal mock communities individually).
Now we will focus on results at species level (for genus level, change to level 5)
End of explanation
for i in [n for n in range(1,27)]:
    display(Markdown('## mock-{0}'.format(i)))
    best = method_by_dataset_a1(mock_results_6, 'mock-{0}'.format(i))
    display(best)
Explanation: Method Optimization
Which method/parameter configuration performed "best" for a given score? We can rank the top-performing configuration by dataset, method, and taxonomic level.
First, the top-performing method/configuration combination by dataset.
End of explanation
for method in mock_results_6['Method'].unique():
    top_params = parameter_comparisons(
        mock_results_6, method, 
        metrics=['Precision', 'Recall', 'F-measure',
                 'Taxon Accuracy Rate', 'Taxon Detection Rate'])
    display(Markdown('## {0}'.format(method)))
    display(top_params[:10])
uniform_6 = mock_results_6[['uniform' in p for p in mock_results_6['Parameters']]]
for method in uniform_6['Method'].unique():
    top_params = parameter_comparisons(
        uniform_6, method, 
        metrics=['Precision', 'Recall', 'F-measure',
                 'Taxon Accuracy Rate', 'Taxon Detection Rate'])
    display(Markdown('## {0}'.format(method)))
    display(top_params[:10])
Explanation: 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).
End of explanation
rank_optimized_method_performance_by_dataset(mock_results,
                                             dataset="Reference",
                                             metric="F-measure",
                                             level_range=range(4,7),
                                             display_fields=["Method",
                                                             "Parameters",
                                                             "Taxon Accuracy Rate",
                                                             "Taxon Detection Rate",
                                                             "Precision",
                                                             "Recall",
                                                             "F-measure"],
                                             paired=True,
                                             parametric=True,
                                             color=None,
                                             color_palette=color_pallette)
rank_optimized_method_performance_by_dataset(mock_results,
                                             dataset="Reference",
                                             metric="Taxon Accuracy Rate",
                                             level_range=range(6,7),
                                             display_fields=["Method",
                                                             "Parameters",
                                                             "Taxon Accuracy Rate",
                                                             "Taxon Detection Rate",
                                                             "Precision",
                                                             "Recall",
                                                             "F-measure"],
                                             paired=True,
                                             parametric=True,
                                             color=None,
                                             color_palette=color_pallette)
Explanation: Optimized method performance
And, finally, which method performed best at each individual taxonomic level for each reference dataset (i.e., for across all fungal and bacterial mock communities combined)?
For this analysis, we rank the top-performing method/parameter combination for each method at family through species levels. Methods are ranked by top F-measure, and the average value for each metric is shown (rather than count best as above). F-measure distributions are plotted for each method, and compared using paired t-tests with FDR-corrected P-values. This cell does not need to be altered, unless if you wish to change the metric used for sorting best methods and for plotting.
End of explanation
<END_TASK>