Community and function¶

This notebook opens the community example, which joins the CAMPI study groups with recorded Unipept taxonomy assignments and one real eggNOG annotation, offline. It shows the three layers FastaLake keeps separate: the taxonomy of accepted peptides, the function of selected proteins, and the peptide to protein evidence network that connects them.

In [1]:

import json
import os
from pathlib import Path

import matplotlib.pyplot as plt
import pandas as pd

# Where the bundled examples wrote their outputs. Point this at your own run to reuse the cells.
DATA = Path(os.environ.get("FASTALAKE_TUTORIAL_DATA", "../../tutorial_data")).resolve()
OKABE_ITO = ["#0072B2", "#E69F00", "#009E73", "#D55E00", "#CC79A7", "#56B4E9", "#F0E442", "#000000"]
plt.rcParams.update({"figure.dpi": 110, "axes.spines.top": False, "axes.spines.right": False, "font.size": 9})
pd.set_option("display.width", 140, "display.max_colwidth", 60)
RUN = DATA / "community"
assert RUN.is_dir(), f"no community run under {DATA}"
print(json.loads((RUN / "VALIDATION.json").read_text()))

import json import os from pathlib import Path import matplotlib.pyplot as plt import pandas as pd # Where the bundled examples wrote their outputs. Point this at your own run to reuse the cells. DATA = Path(os.environ.get("FASTALAKE_TUTORIAL_DATA", "../../tutorial_data")).resolve() OKABE_ITO = ["#0072B2", "#E69F00", "#009E73", "#D55E00", "#CC79A7", "#56B4E9", "#F0E442", "#000000"] plt.rcParams.update({"figure.dpi": 110, "axes.spines.top": False, "axes.spines.right": False, "font.size": 9}) pd.set_option("display.width", 140, "display.max_colwidth", 60) RUN = DATA / "community" assert RUN.is_dir(), f"no community run under {DATA}" print(json.loads((RUN / "VALIDATION.json").read_text()))

{'status': 'PASS', 'samples': 2, 'network_access': False, 'taxonomy': 'Recorded native Unipept output; reference-dependent LCA', 'function': 'One real sequence-matched enolase annotation; partial coverage', 'ckg_renderer': 'Separate optional environment; not executed by this example', 'counts_and_signal_conserved': True}

Peptide taxonomy¶

taxonomy/peptide_taxonomy.tsv carries Unipept's lowest common ancestor per accepted peptide, with the summed LFQ intensity behind it. The rank at which a peptide resolves is itself a result: many peptides stop at class or phylum.

In [2]:

tax = pd.read_csv(RUN / "taxonomy/peptide_taxonomy.tsv", sep="\t")
print(tax["status"].value_counts().to_dict())
rank_order = ["superkingdom", "phylum", "class", "order", "family", "genus", "species"]
ranks = tax[tax["status"] == "assigned"].groupby(["sample", "taxon_rank"]).size().unstack(fill_value=0)
ranks[[r for r in rank_order if r in ranks.columns]]

tax = pd.read_csv(RUN / "taxonomy/peptide_taxonomy.tsv", sep="\t") print(tax["status"].value_counts().to_dict()) rank_order = ["superkingdom", "phylum", "class", "order", "family", "genus", "species"] ranks = tax[tax["status"] == "assigned"].groupby(["sample", "taxon_rank"]).size().unstack(fill_value=0) ranks[[r for r in rank_order if r in ranks.columns]]

{'assigned': 751, 'unresolved': 218, 'lookup_cutoff': 5}

Out[2]:

taxon_rank	phylum	class	order	family	genus
sample
S01	21	32	59	69	112
S02	21	22	44	45	94

In [3]:

genus = tax[(tax["status"] == "assigned") & (tax["taxon_rank"] == "genus")]
top = genus.groupby(["taxon_name", "sample"])["lfq_sum_intensity"].sum().unstack(fill_value=0)
top = top.loc[top.sum(axis=1).sort_values(ascending=False).index[:10]]
ax = top.plot.barh(figsize=(6, 3.4), color=OKABE_ITO[:top.shape[1]], width=0.8)
ax.invert_yaxis(); ax.set_xlabel("summed LFQ intensity of genus-resolved peptides"); ax.set_ylabel(""); ax.legend(frameon=False)
plt.tight_layout()

genus = tax[(tax["status"] == "assigned") & (tax["taxon_rank"] == "genus")] top = genus.groupby(["taxon_name", "sample"])["lfq_sum_intensity"].sum().unstack(fill_value=0) top = top.loc[top.sum(axis=1).sort_values(ascending=False).index[:10]] ax = top.plot.barh(figsize=(6, 3.4), color=OKABE_ITO[:top.shape[1]], width=0.8) ax.invert_yaxis(); ax.set_xlabel("summed LFQ intensity of genus-resolved peptides"); ax.set_ylabel(""); ax.legend(frameon=False) plt.tight_layout()

holobiont_balance.tsv splits the same signal into host, bacteria, archaea, eukaryote and unassigned fractions, by peptide count and by intensity.

In [4]:

balance = pd.read_csv(RUN / "taxonomy/holobiont_balance.tsv", sep="\t")
balance.pivot(index="category", columns="sample", values=["peptides", "signal_fraction"]).round(3)

balance = pd.read_csv(RUN / "taxonomy/holobiont_balance.tsv", sep="\t") balance.pivot(index="category", columns="sample", values=["peptides", "signal_fraction"]).round(3)

Out[4]:

	peptides		signal_fraction
sample	S01	S02	S01	S02
category
Archaea	0.0	0.0	0.000	0.000
Bacteria	420.0	327.0	0.740	0.806
Host	0.0	0.0	0.000	0.000
Other assigned taxa	3.0	1.0	0.002	0.001
Unresolved	125.0	98.0	0.258	0.193

Function of selected proteins¶

functions/protein_annotations.tsv has one row per selected protein with its annotation status. This offline example carries exactly one real sequence-matched eggNOG annotation (an enolase); every other protein is explicitly not_annotated. Partial coverage is stated, never filled in.

In [5]:

ann = pd.read_csv(RUN / "functions/protein_annotations.tsv", sep="\t")
print(ann["annotation_status"].value_counts().to_dict())
ann[ann["annotation_status"] == "assigned"][["protein_id", "match_rule", "KEGG_ko", "EC", "PFAMs"]]

ann = pd.read_csv(RUN / "functions/protein_annotations.tsv", sep="\t") print(ann["annotation_status"].value_counts().to_dict()) ann[ann["annotation_status"] == "assigned"][["protein_id", "match_rule", "KEGG_ko", "EC", "PFAMs"]]

{'not_annotated': 447, 'assigned': 1}

Out[5]:

	protein_id	match_rule	KEGG_ko	EC	PFAMs
112	CAMPI_427a9c2476c475ddd52e2c34117aff61006a9c1cc1efc2ed5e...	exact_sequence_alias	ko:K01689	4.2.1.11	Enolase_C,Enolase_N

The evidence network¶

network/nodes.tsv and network/edges.tsv describe which peptide was reported against which protein, in which acquisitions, and which protein is the assigned representative. A peptide with reported degree above one is the shared evidence that keeps a group ambiguous.

In [6]:

nodes = pd.read_csv(RUN / "network/nodes.tsv", sep="\t")
edges = pd.read_csv(RUN / "network/edges.tsv", sep="\t")
print(nodes["node_type"].value_counts().to_dict(), "| edges:", len(edges))
peptides = nodes[nodes["node_type"] == "peptide"]
print("peptides by reported degree:", peptides["reported_degree"].value_counts().sort_index().to_dict())
print("peptides observed in both acquisitions:", int((peptides["observed_acquisitions"] == 2).sum()))
peptides[peptides["reported_degree"] > 1][["label", "reported_degree", "observed_acquisitions", "assigned_group"]]

nodes = pd.read_csv(RUN / "network/nodes.tsv", sep="\t") edges = pd.read_csv(RUN / "network/edges.tsv", sep="\t") print(nodes["node_type"].value_counts().to_dict(), "| edges:", len(edges)) peptides = nodes[nodes["node_type"] == "peptide"] print("peptides by reported degree:", peptides["reported_degree"].value_counts().sort_index().to_dict()) print("peptides observed in both acquisitions:", int((peptides["observed_acquisitions"] == 2).sum())) peptides[peptides["reported_degree"] > 1][["label", "reported_degree", "observed_acquisitions", "assigned_group"]]

{'peptide': 819, 'reported_protein': 448} | edges: 826
peptides by reported degree: {1: 814, 2: 3, 3: 2}
peptides observed in both acquisitions: 155

Out[6]:

	label	reported_degree	observed_acquisitions	assigned_group
181	FQHPVAGTYK	3	2	CAMPI_8b9961cc9da37de018df0b80a1b363618ea56e109ff83e0ff4...
182	FQHPVAGTYKK	2	1	CAMPI_8b9961cc9da37de018df0b80a1b363618ea56e109ff83e0ff4...
277	HYAHVDCPGHADYVK	3	2	CAMPI_c441c84224d4be45b80629c2a8a57b3eefbb8b205bc299e5c2...
301	KGLADTALK	2	1	CAMPI_7cb1e4d122178b054c32ee50ef164d0d60b9618226ddad845b...
342	KYFSVASGGGTGR	2	1	CAMPI_8b9961cc9da37de018df0b80a1b363618ea56e109ff83e0ff4...