"""Helpers to convert external data to a MICOM taxonomy."""
from micom.qiime_formats import (
load_qiime_feature_table,
load_qiime_taxonomy,
)
import pandas as pd
[docs]
RANKS = ["kingdom", "phylum", "class", "order", "family", "genus", "species", "strain"]
[docs]
def build_from_qiime(
abundance: pd.DataFrame,
taxonomy: pd.Series,
collapse_on: str = "genus",
trim_rank_prefix: bool = False,
) -> pd.DataFrame:
"""Build the specification for the community models."""
if trim_rank_prefix:
taxa = taxonomy.str.replace("[\\w_]+__|\\[|\\]", "", regex=True)
else:
taxa = taxonomy
taxa = taxa.str.split(";\\s*", expand=True).replace("", None)
taxa.columns = RANKS[0 : taxa.shape[1]]
taxa["taxid"] = taxonomy.index
taxa.index == taxa.taxid
if isinstance(collapse_on, str):
collapse_on = [collapse_on]
ranks = [r for r in collapse_on if r in taxa.columns]
taxa["mapping_ranks"] = taxa[ranks].apply(
lambda s: "|".join(s.astype("str")), axis=1
)
abundance = (
abundance.collapse(
lambda id_, x: taxa.loc[id_, "mapping_ranks"],
axis="observation",
norm=False,
)
.to_dataframe(dense=True)
.T
)
abundance["sample_id"] = abundance.index
abundance = abundance.melt(
id_vars="sample_id", var_name="mapping_ranks", value_name="abundance"
)
abundance = pd.merge(
abundance[abundance.abundance > 0.0],
taxa[ranks + ["mapping_ranks"]].drop_duplicates(),
on="mapping_ranks",
)
abundance["id"] = abundance["mapping_ranks"].replace(
r"[^A-Za-z0-9_]+", "_", regex=True
)
del abundance["mapping_ranks"]
abundance.dropna(subset=ranks, inplace=True)
depth = abundance.groupby("sample_id").abundance.sum()
abundance["relative"] = abundance.abundance / depth[abundance.sample_id].values
return abundance
[docs]
def qiime_to_micom(feature_table, taxonomy, collapse_on="genus"):
"""Load a micom taxonomy from Qiime 2 data.
Parameters
----------
feature_table : str
Path to a Qiime 2 FeatureTable artifact.
taxonomy : str
Path to a Qiime 2 FeatureData[Taxonomy] artifact.
collapse_on : str or List[str]
The taxa ranks to collapse on. This will dictate how strict the database
matching will be as well.
Returns
-------
pd.DataFrame
A micom taxonomy containing abundances and taxonomy calls in long
format.
"""
table = load_qiime_feature_table(feature_table)
taxonomy = load_qiime_taxonomy(taxonomy)
return build_from_qiime(table, taxonomy, collapse_on)
[docs]
def rank_prefixes(manifest: pd.DataFrame) -> pd.Series:
"""Get the used prefixes for taxonomic ranks.
Arguments
---------
manifest : pandas.DataFrame
A model database manifest.
Returns
-------
pandas.Series
The detected prefix for each taxonomic rank in the manifest.
"""
ranks = [c for c in manifest.columns if c.lower() in RANKS]
prefixes = pd.Series(
{r: manifest[r].str.extract(r"^([a-z]__)").iloc[0, 0] for r in ranks}
)
return prefixes
[docs]
def unify_rank_prefixes(taxonomy: pd.DataFrame, manifest: pd.DataFrame) -> pd.DataFrame:
"""Handle taxonomic rank prefixes in the taxonomy or database manifest.
Arguments
---------
taxonomy : pandas.DataFrame
A taxonomy table.
manifest : pandas.DataFrame
A database manifest.
Returns
-------
tuple of pandas.DataFrame
The taxonomy with adjusted taxa names consistent with the database.
"""
tax_prefixes = rank_prefixes(taxonomy)
db_prefixes = rank_prefixes(manifest)
ranks = tax_prefixes.index[tax_prefixes.index.isin(db_prefixes.index)]
if all(tax_prefixes[ranks] == db_prefixes[ranks]):
return taxonomy
taxonomy = taxonomy.copy()
ranks = [c for c in taxonomy.columns if c.lower() in RANKS]
if db_prefixes.isna().all():
for r in ranks:
taxonomy[r] = taxonomy[r].str.replace(r"^[a-z]__", "", regex=True)
else:
for r in ranks:
taxonomy[r] = db_prefixes[r] + taxonomy[r]
return taxonomy