<jupyter_start><jupyter_text>Figure 1A - Full EMP VisualizationThis will walkthrough the procedure of generating the EMPire plot of the full EMP dataset.<jupyter_code># Clear the output directory and create a fresh one
!rm -rf fig1/output
!mkdir -p fig1/output<jupyter_output><empty_output><jupyter_text>Downloading Required DataFirst, we will download the EMP data via FTP. Note that we do not include these files in the repository because they are fairly large.<jupyter_code># EMP table
!wget -O fig1/input/EMP_tbl.biom ftp://ftp.microbio.me/emp/release1/otu_tables/deblur/emp_deblur_90bp.release1.biom
# EMP tree
!wget -O fig1/input/EMP_tree.tre ftp://ftp.microbio.me/emp/release1/otu_info/deblur/emp90.5000_1000_rxbl_placement_pruned75.tog.tre
# EMP mapping file
!wget -O fig1/input/EMP_metadata.tsv ftp://ftp.microbio.me/emp/release1/mapping_files/emp_qiime_mapping_release1.tsv<jupyter_output>--2020-10-01 10:38:18--  ftp://ftp.microbio.me/emp/release1/mapping_files/emp_qiime_mapping_release1.tsv
           => ‘fig1/input/EMP_metadata.tsv’
Resolving ftp.microbio.me (ftp.microbio.me)... 192.168.3.11
Connecting to ftp.microbio.me (ftp.microbio.me)|192.168.3.11|:21... connected.
Logging in as anonymous ... Logged in!
==> SYST ... done.    ==> PWD ... done.
==> TYPE I ... done.  ==> CWD (1) /emp/release1/mapping_files ... done.
==> SIZE emp_qiime_mapping_release1.tsv ... 21876184
==> PASV ... done.    ==> RETR emp_qiime_mapping_release1.tsv ... done.
Length: 21876184 (21M) (unauthoritative)

emp_qiime_mapping_r 100%[===================>]  20.86M  16.9MB/s    in 1.2s    

2020-10-01 10:38:20 (16.9 MB/s) - ‘fig1/input/EMP_metadata.tsv’ saved [21876184]<jupyter_text>Next, we'll import the table and tree files as `qza` Artifacts so that they can be used with EMPress.<jupyter_code>!qiime tools import \
    --type FeatureTable[Frequency] \
    --input-path fig1/input/EMP_tbl.biom \
    --output-path fig1/output/EMP_tbl.qza
!qiime tools import \
    --type Phylogeny[Rooted] \
    --input-path fig1/input/EMP_tree.tre \
    --output-path fig1/output/EMP_tree.qza<jupyter_output>[32mImported fig1/input/EMP_tree.tre as NewickDirectoryFormat to fig1/output/EMP_tree.qza[0m<jupyter_text>Simplify Phylum-Level AssignmentsFor ease of visualization, we will use each microbial feature's phylum level assignment for subsequent taxonomic representation. These assignments are stored within the original `biom` table. We will first extract this file and replace missing values with "Unspecified."<jupyter_code>tbl_biom = load_table("fig1/input/EMP_tbl.biom")
taxonomy_df = tbl_biom.metadata_to_dataframe(axis="observation")
#taxonomy_df.columns = ["kingdom", "phylum", "class", "order", "family", "genus", "species"]
taxonomy_df.index.name = "featureid"
taxonomy_df = taxonomy_df.fillna("Unspecified")
print(taxonomy_df.shape)
taxonomy_df.head()<jupyter_output>(317314, 7)<jupyter_text>We will save this taxonomy file as its own feature metadata file independent of the feature table.<jupyter_code>Metadata(taxonomy_df).save("fig1/output/EMP_taxonomy.tsv")<jupyter_output><empty_output><jupyter_text>We want to keep only the top 5 most abundant phlyum (taxonomy_1) assignments.<jupyter_code>taxonomy_df["taxonomy_1"].value_counts().head(6) # top 5 + Unspecified
top_5_phyla_counts = taxonomy_df["taxonomy_1"].value_counts()[:6]
top_5_phyla = top_5_phyla_counts.index.tolist()
bot_phyla = taxonomy_df["taxonomy_1"].value_counts()[6:].index.tolist()
top_5_phyla.remove("Unspecified")
bot_phyla.append("Unspecified")
phylum_assignments = ["Other/Unspecified" if not x in top_5_phyla else x for x in taxonomy_df["taxonomy_1"]]
top_5_phyla_df = pd.DataFrame.from_dict(dict(zip(taxonomy_df.index.tolist(), phylum_assignments)), orient="index")
top_5_phyla_df.columns = ["phylum_assignment"]
top_5_phyla_df.index.name = "featureid"
top_5_phyla_df.head()<jupyter_output><empty_output><jupyter_text>We will save this file as a `tsv` file that QIIME2 can read.<jupyter_code>Metadata(top_5_phyla_df).save("fig1/output/EMP_top_5_phyla.tsv")<jupyter_output><empty_output><jupyter_text>Filtering Feature TableThe original EMP data includes samples from both sterile water blanks and mock communities. We will remove samples of both these types.<jupyter_code>!qiime feature-table filter-samples \
    --i-table fig1/output/EMP_tbl.qza \
    --m-metadata-file fig1/input/EMP_metadata.tsv \
    --p-where "[empo_2] NOT IN ('Positive', 'Negative')" \
    --o-filtered-table fig1/output/EMP_tbl_filt.qza<jupyter_output>[32mSaved FeatureTable[Frequency] to: fig1/output/EMP_tbl_filt.qza[0m<jupyter_text>Additionally, we will filter the feature table by the phylogeny to make sure all features present in the feature table are present in the phylogeny as well.<jupyter_code>!qiime phylogeny filter-table \
    --i-table fig1/output/EMP_tbl_filt.qza \
    --i-tree fig1/output/EMP_tree.qza \
    --o-filtered-table fig1/output/EMP_tbl_filt.qza<jupyter_output>[32mSaved FeatureTable[Frequency] to: fig1/output/EMP_tbl_filt.qza[0m<jupyter_text>Filtering MetadataThe EMP dataset is fairly large so the metadata takes up a good amount of storage. We will pare down the provided metadata file a bit so that the final Visualization is small enough to upload to GitHub.<jupyter_code>metadata = pd.read_csv("fig1/input/EMP_metadata.tsv", sep="\t", index_col=0)
metadata.head()<jupyter_output><empty_output><jupyter_text>We are going to keep only a few metadata categories that we think are useful for interactive visualization.<jupyter_code>cols_to_keep = [
    "empo_1", "empo_2", "empo_3", "study_id", "sample_scientific_name",
    "host_common_name", "env_material", "country"
]
metadata_filt = metadata[cols_to_keep]
print(metadata_filt.shape)
Metadata(metadata_filt).save("fig1/input/EMP_metadata_filt.tsv")<jupyter_output>(27751, 8)<jupyter_text>Running EMPressFinally, we will generate an EMPire plot combining an interactive tree visualization with a PCoA. The distance matrix and PCoA are very computationally intensive to calculate so we have included the pre-computed PCoA in `fig1/input`. The code used to generate the unweighted UniFrac distance matrix and PCoA are shown below (not executed). ```bashqiime diversity beta-phylogenetic \    --i-table fig1/ouptut/EMP_tbl_filt.qza \    --i-phylogeny fig1/output/EMP_tree.qza \    --p-metric unweighted_unifrac \    --p-n-jobs 32 \    --o-distance-matrix fig1/output/EMP_u_unifrac_dm.qza    qiime diversity pcoa \    --i-distance-matrix fig1/output/EMP_u_unifrac_dm.qza \    --p-number-of-dimensions 3 \    --o-pcoa fig1/input/EMP_u_unifrac_pcoa.qza```<jupyter_code>!qiime empress community-plot \
    --i-tree fig1/output/EMP_tree.qza \
    --i-feature-table fig1/output/EMP_tbl_filt.qza \
    --i-pcoa fig1/input/EMP_u_unifrac_pcoa.qza \
    --m-sample-metadata-file fig1/input/EMP_metadata_filt.tsv \
    --m-feature-metadata-file fig1/output/EMP_top_5_phyla.tsv \
    --m-feature-metadata-file fig1/output/EMP_taxonomy.tsv \
    --o-visualization fig1/output/EMP_empire.qzv<jupyter_output>[32mSaved Visualization to: fig1/output/EMP_empire.qzv[0m<jupyter_text>Figure 1BNow we will recreate the subtree with the pH metadata. First, we will import the sample metadata to see the overall distribution of pH values.<jupyter_code>ph_vals = metadata["ph"].dropna()
plt.boxplot(ph_vals)
plt.show()<jupyter_output><empty_output><jupyter_text>To remove outliers, we are going to retain only the samples in the middle 90% of the range of pH values.<jupyter_code>q1, q2 = [ph_vals.quantile(x) for x in [0.05, 0.95]]
print(q1, q2)
!qiime feature-table filter-samples \
    --i-table fig1/output/EMP_tbl_filt.qza \
    --m-metadata-file fig1/input/EMP_metadata.tsv \
    --p-where "[ph] >= 4.7 AND [ph] <= 9" \
    --o-filtered-table fig1/output/ph_EMP_tbl_filt.qza<jupyter_output>[32mSaved FeatureTable[Frequency] to: fig1/output/ph_EMP_tbl_filt.qza[0m<jupyter_text>For each feature, we will calculate:1. The number of samples in which it appears2. The arithmetic mean pH value of these samples3. The geometric mean pH value of these samples<jupyter_code>ph_tbl_df = Artifact.load("fig1/output/ph_EMP_tbl_filt.qza").view(pd.DataFrame)
features = ph_tbl_df.columns.tolist()
feature_ph_dict = dict()

for feature in features:
    feature_ph_dict[feature] = dict()
    samps = ph_tbl_df.loc[:, feature]
    samps = samps[samps > 0].index.tolist()
    ph_vals = metadata.loc[samps]["ph"].tolist()
    feature_ph_dict[feature]["sample_presence"] = len(samps)
    feature_ph_dict[feature]["geometric_mean"] = ss.mstats.gmean(ph_vals)
    feature_ph_dict[feature]["arithmetic_mean"] = np.mean(ph_vals)
    
ph_mean_df = pd.DataFrame.from_dict(feature_ph_dict, orient="index")
ph_mean_df.index.name = "featureid"
ph_mean_df.head()
# log10 samples to tighten scale for colorbar
ph_mean_df["log10_sample_presence"] = ph_mean_df["sample_presence"].apply(np.log10)
Metadata(ph_mean_df).save("fig1/output/ph_EMP_feature_metadata.tsv")
!qiime empress community-plot \
    --i-tree fig1/output/EMP_tree.qza \
    --i-feature-table fig1/output/ph_EMP_tbl_filt.qza \
    --m-sample-metadata-file fig1/input/EMP_metadata_filt.tsv \
    --m-feature-metadata-file fig1/output/EMP_top_5_phyla.tsv \
    --m-feature-metadata-file fig1/output/ph_EMP_feature_metadata.tsv \
    --m-feature-metadata-file fig1/output/EMP_taxonomy.tsv \
    --o-visualization fig1/output/ph_EMP_empress.qzv<jupyter_output>[32mSaved Visualization to: fig1/output/ph_EMP_empress.qzv[0m<jupyter_text>Figure 1CFinally, we will plot the distributions of pH for the top 5 most abundant phyla in the EMP features.<jupyter_code>ph_taxonomy = taxonomy_df.join(ph_mean_df)
ph_taxonomy.head()<jupyter_output><empty_output><jupyter_text>We will first perform a similar procedure as above where we mask phyla outside the top 5 most abundant as "Other/Unspecified". Then we will assign the order of the top 5 phyla in increasing median mean pH.<jupyter_code>other_phyla = [x for x in ph_taxonomy["taxonomy_1"].unique() if x not in top_5_phyla]
ph_phyla_mean_df = ph_taxonomy.copy()[["taxonomy_1", "arithmetic_mean"]]
ph_phyla_mean_df["taxonomy_1"] = ph_phyla_mean_df["taxonomy_1"].replace({x: "Other/Unspecified" for x in other_phyla})
ph_phyla_mean_df = ph_phyla_mean_df.dropna()
order = ph_phyla_mean_df.groupby("taxonomy_1").median().sort_values(by="arithmetic_mean", ascending=True).index.tolist()
order.remove("Other/Unspecified") # want Other/Unspecified to appear last
order.append("Other/Unspecified")
fig, axs = plt.subplots(6, 1, figsize=(6, 8), sharex=True)
plt.subplots_adjust(hspace=0.1)
order_dict = dict(zip(order, np.arange(len(order))))

ph_phyla_mean_df["rank"] = ph_phyla_mean_df["taxonomy_1"].map(order_dict)
ph_phyla_mean_df = ph_phyla_mean_df.sort_values(by="rank", ascending=True)

palette = dict(zip(
    [       
        "Other/Unspecified",
        "p__Bacteroidetes",
        "p__Chloroflexi",
        "p__Firmicutes",
        "p__Planctomycetes",
        "p__Proteobacteria",
    ],
    sns.color_palette("Dark2", 6))
)

for i, ax in enumerate(axs.flatten()):
    p = order[i]
    vals = ph_phyla_mean_df[ph_phyla_mean_df["taxonomy_1"] == p]["arithmetic_mean"]
    _g = ax.hist(
        vals,
        bins=50,
        color=palette[p],
    )
    ax.axvline(x = np.median(vals), color="black", linestyle="--")
    if i != 5: ax.xaxis.set_visible(False)
        
    _name = p
    if _name.startswith("p__"): _name = _name[3:]
    ax.text(
        x=0.05,
        y=0.9,
        s=f"{_name}",
        color=palette[p],
        ha="left",
        va="top",
        fontweight="bold",
        transform=ax.transAxes,
        fontsize=14,
    )
    
    _y1, _y2 = ax.get_ylim()    
    ax.set_ylim([_y1, _y2*1.5])
    
    ax.set_yticks([0, _g[0].max()])
    ax.set_yticklabels(map(int, [0, _g[0].max()]), fontsize=14)
    ax.tick_params("x", labelsize=14)

fig.add_subplot(111, frameon=False)
# hide tick and tick label of the big axes
plt.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
plt.grid(False)
plt.ylabel("Number of Features", fontsize=20, labelpad=20)
plt.xlabel("Mean pH", fontsize=20)

plt.savefig("fig1/output/mean_ph_hist.svg", dpi=300, bbox_inches="tight")

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12647898816809466, 'fraction_numerical': 0.025785393716850267, 'mean_word_length': 4.794799054373523, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 71, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '423983', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4769, 'n_tokens_neox': 4546, 'n_words': 906, 'path': 'notebooks/fig1.ipynb', 'repo_name': 'kwcantrell/empress-analyses', 'size': '12257'}