<jupyter_start><jupyter_text>Greengenes OTU Clustering BenchmarkThis notebook tests the effects of clustering at differnt OTU thresholds on taxonomic information and classification accuracy, using the Greengenes 16S rRNA gene reference database.<jupyter_code>import os
import qiime2 as q2
from qiime2.plugins import rescript
from glob import glob
import timeit
import tempfile
from glob import glob
import seaborn as sns
from matplotlib import pyplot as plt

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Define inputsModify the filepaths here to your local filepaths to run this notebook.The filepaths below merely point to the sequence and taxonomy files for the greengenes 13_8 release, clustered at different OTU % similarities.Greengenes files can be downloaded and unzipped from here:ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz**DO NOT MODIFY ANY OTHER CELLS IN THIS NOTEBOOK IF YOU ARE ATTEMPTING TO REPLICATE THESE RESULTS**<jupyter_code># N jobs to use
n_jobs = 4

# location of reference database files
dbdir = '../../../ref_dbs/gg_13_8_otus/'
dbs = {
    'gg_{0}'.format(n): {
        'seqs': q2.Artifact.import_data(
            'FeatureData[Sequence]', dbdir + 'rep_set/{0}_otus.fasta'.format(n)),
        'taxa': q2.Artifact.import_data(
            'FeatureData[Taxonomy]', dbdir + 'taxonomy/{0}_otu_taxonomy.txt'.format(n),
            view_type='HeaderlessTSVTaxonomyFormat')}
    for n in ['64', '79', '88', '91', '94', '97', '99']}<jupyter_output><empty_output><jupyter_text>Perform CV ClassificationBoth k-fold and "perfect classifiers" (no CV, leaky data, best possible performance)<jupyter_code># Cross-validation
times_cv = {}
for n, v in dbs.items():
    s, t = v['seqs'], v['taxa']
    print(n)
    exp_taxa_fp = '../data/expected/{0}_expected_taxonomy.qza'.format(n)
    if not os.path.exists(exp_taxa_fp):
        start = timeit.default_timer()
        exp, obs, = rescript.actions.evaluate_cross_validate(s, t, n_jobs=n_jobs)
        times_cv[n] = timeit.default_timer() - start
        exp.save(exp_taxa_fp)
        obs.save('../data/observed/{0}_observed_taxonomy.qza'.format(n))
# "Perfect" classifier, no cross-validation
times_nocv = {}
for n, v in dbs.items():
    s, t = v['seqs'], v['taxa']
    print(n)
    obs_taxa_fp = '../data/observed/{0}_observed_taxonomy_noCV.qza'.format(n)
    exp_taxa_fp = '../data/expected/{0}_expected_taxonomy_noCV.qza'.format(n)
    if not os.path.exists(exp_taxa_fp):
        start = timeit.default_timer()
        _, _, obs, = rescript.actions.evaluate_fit_classifier(s, t, n_jobs=n_jobs)
        times_nocv[n] = timeit.default_timer() - start
        obs.save(obs_taxa_fp)
        t.save(exp_taxa_fp)

exp_taxonomies = {os.path.splitext(os.path.basename(fp))[0].replace('_expected_taxonomy', ''): 
                  q2.Artifact.load(fp) for fp in glob('../data/expected/gg_*')}
obs_taxonomies = {os.path.splitext(os.path.basename(fp))[0].replace('_observed_taxonomy', ''): 
                  q2.Artifact.load(fp) for fp in glob('../data/observed/gg_*')}<jupyter_output><empty_output><jupyter_text>Evaluate Classification Accuracy<jupyter_code>class_eval, = rescript.actions.evaluate_classifications(
    expected_taxonomies = [exp_taxonomies[i] for i in sorted(exp_taxonomies.keys())],
    observed_taxonomies = [obs_taxonomies[i] for i in sorted(obs_taxonomies.keys())],
    labels = sorted(obs_taxonomies.keys()))
class_eval.save('../data/eval/gg_otu_cluster_classification_eval.qzv')<jupyter_output><empty_output><jupyter_text>Evaluate Taxonomic Information<jupyter_code>tax_eval, = rescript.actions.evaluate_taxonomy(
    taxonomies = [dbs[i]['taxa'] for i in sorted(dbs.keys())],
    labels = sorted(dbs.keys()))
tax_eval.save('../data/eval/gg_otu_cluster_taxonomic_eval.qzv')
tax_eval, = rescript.actions.evaluate_taxonomy(
    taxonomies = [dbs[i]['taxa'] for i in sorted(dbs.keys())],
    labels = sorted(dbs.keys()),
    rank_handle_regex="^[dkpcofgs]__")
tax_eval.save('../data/eval/gg_otu_cluster_taxonomic_eval_dropemptyranks.qzv')<jupyter_output><empty_output><jupyter_text>Evaluate Sequence Information<jupyter_code>seqs_eval, = rescript.actions.evaluate_seqs(
    sequences = [dbs[i]['seqs'] for i in sorted(dbs.keys())],
    labels = sorted(dbs.keys()),
    palette = 'cividis')
seqs_eval.save('../data/eval/gg_otu_cluster_seqs_eval.qzv')<jupyter_output><empty_output><jupyter_text>Make a pretty figureThe interactive figs are great, but we need to generate a static figure for publication :(<jupyter_code>fig_data = dict()
with tempfile.TemporaryDirectory() as temp_dir_name:
    for f in glob('../data/eval/*.qzv'):
        name = f.split('/')[-1].split('.qzv')[0]
        f = q2.Visualization.load(f)
        f.export_data(temp_dir_name)
        fig_data[name] = q2.Metadata.load(temp_dir_name + '/data.tsv').to_dataframe()

plot_keys = ['gg_otu_cluster_taxonomic_eval',
             'gg_otu_cluster_taxonomic_eval',
             'gg_otu_cluster_taxonomic_eval_dropemptyranks',
             'gg_otu_cluster_classification_eval',
             'gg_otu_cluster_classification_eval']
signals = ['Unique Labels',
           'Taxonomic Entropy',
           'Proportion of Features Terminating at Depth',
           'F-Measure',
           'F-Measure']
labels = list('ABCDEF')

fig, axes = plt.subplots(2, 3, figsize=(12, 6), sharex=False)
axes = axes.flatten()
axes[-1].remove()
for ax, key, signal, label in zip(axes, plot_keys, signals, labels):
    ax.text(-0.1, 1.1, label, transform=ax.transAxes,
            fontsize=16, fontweight='bold', va='top', ha='right')
    dat = fig_data[key].copy()
    dat['Level'] = dat['Level'].replace({
        1: 'D', 2: 'P', 3: 'C', 4: 'O',
        5: 'F', 6: 'G', 7: 'S'})
    if ax == axes[3]:
        dat = dat[dat['Dataset'].isin(['gg_64_noCV', 'gg_79_noCV', 'gg_88_noCV',
                                       'gg_91_noCV', 'gg_94_noCV', 'gg_97_noCV',
                                       'gg_99_noCV'])]
    if ax == axes[4]:
        dat = dat[dat['Dataset'].isin(['gg_64', 'gg_79', 'gg_88', 'gg_91',
                                       'gg_94', 'gg_97', 'gg_99'])]
        legend = 'brief'
    else:
        legend = False
    g = sns.lineplot(x='Level', y=signal, hue='Dataset', data=dat, ax=ax,
                     legend=legend, palette='cividis', sort=False)
    if ax == axes[2]:
        ax.set_ylabel('Terminal Labels')
    if ax == axes[4]:
        ax.legend(bbox_to_anchor=(1.1, 1.05))
    else:
        ax.set_xlabel('')
fig.savefig('../data/plots/greengenes_evaluation.pdf', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1445619335347432, 'fraction_numerical': 0.015861027190332326, 'mean_word_length': 3.804789550072569, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '566014', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2461, 'n_tokens_neox': 2334, 'n_words': 450, 'path': 'notebooks/gg-cluster-test.ipynb', 'repo_name': 'bokulich-lab/db-benchmarks-2020', 'size': '6620'}