Matched DNA and RNA¶
This notebook opens the multi-omics example: the real CAMPI spectra with synthetic metagenome and metatranscriptome values. It shows the ledger that explains why every molecular candidate was included or excluded, what the additions changed in the search, and how the two searches are compared. The TPM values are artificial test inputs; the mechanics are the real ones.
In [1]:
Copied!
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 / "multiomics"
assert RUN.is_dir(), f"no multi-omics run under {DATA}"
v = json.loads((RUN / "VALIDATION.json").read_text())
print(v["status"]); pd.DataFrame(v["samples"])
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 / "multiomics"
assert RUN.is_dir(), f"no multi-omics run under {DATA}"
v = json.loads((RUN / "VALIDATION.json").read_text())
print(v["status"]); pd.DataFrame(v["samples"])
PASS
Out[1]:
| sample | baseline_sequences | additions | final_sequences | fasta_sha256 | |
|---|---|---|---|---|---|
| 0 | S01 | 348 | 3 | 351 | ae83277ca1dd562d276247df9f83f1216bbe41bc899054ab8158be09... |
| 1 | S02 | 287 | 3 | 290 | c2c073188bd918626c7cbfc7118994e004427810d46bc1bc9cd7a47e... |
The molecular inputs¶
Per specimen: a protein FASTA, a TPM table with metaG_tpm and metaT_tpm (NA where a layer was not measured), and a provenance record.
In [2]:
Copied!
tpm = pd.read_csv(RUN / "inputs/S01/tpm.tsv", sep="\t")
print(json.loads((RUN / "inputs/S01/provenance.json").read_text()))
tpm.head()
tpm = pd.read_csv(RUN / "inputs/S01/tpm.tsv", sep="\t")
print(json.loads((RUN / "inputs/S01/provenance.json").read_text()))
tpm.head()
{'kind': 'synthetic molecular values and annotations over real CAMPI proteins', 'spectra': 'measured bundled CAMPI acquisition window', 'biological_performance_claim': False}
Out[2]:
| prodigal_protein | metaG_tpm | metaT_tpm | |
|---|---|---|---|
| 0 | g0 | 100 | NaN |
| 1 | g1 | 0 | 100.0 |
| 2 | g2 | 0 | NaN |
| 3 | g3 | 0 | NaN |
| 4 | g4 | 0 | NaN |
The gene evidence ledger¶
workflow/samples/<sample>/molecular/gene_evidence.tsv has one row per candidate gene: its TPM values, which thresholds it passed, whether its sequence was already a de novo target, and whether it was included. Every de novo target stays; molecular candidates are added, never substituted.
In [3]:
Copied!
ledger = pd.read_csv(RUN / "workflow/samples/S01/molecular/gene_evidence.tsv", sep="\t")
print(len(ledger), "candidate genes")
pd.crosstab([ledger["dna_pass"], ledger["rna_pass"]], [ledger["in_baseline"], ledger["included"]],
rownames=["dna_pass", "rna_pass"], colnames=["in_baseline", "included"])
ledger = pd.read_csv(RUN / "workflow/samples/S01/molecular/gene_evidence.tsv", sep="\t")
print(len(ledger), "candidate genes")
pd.crosstab([ledger["dna_pass"], ledger["rna_pass"]], [ledger["in_baseline"], ledger["included"]],
rownames=["dna_pass", "rna_pass"], colnames=["in_baseline", "included"])
351 candidate genes
Out[3]:
| in_baseline | 0 | 1 | |
|---|---|---|---|
| included | 1 | 1 | |
| dna_pass | rna_pass | ||
| 0 | 0 | 1 | 348 |
| 1 | 1 | 0 | |
| 1 | 0 | 1 | 0 |
In [4]:
Copied!
ledger[(ledger["included"] == 1) & (ledger["in_baseline"] == 0)][["gene_id", "metaG_tpm", "metaT_tpm", "dna_pass", "rna_pass", "final_accession"]]
ledger[(ledger["included"] == 1) & (ledger["in_baseline"] == 0)][["gene_id", "metaG_tpm", "metaT_tpm", "dna_pass", "rna_pass", "final_accession"]]
Out[4]:
| gene_id | metaG_tpm | metaT_tpm | dna_pass | rna_pass | final_accession | |
|---|---|---|---|---|---|---|
| 0 | g0 | 100 | NaN | 1 | 0 | FMO_SHA256_001f1586f7defaa3805fcb4cf73f6597963087434516e... |
| 1 | g1 | 0 | 100.0 | 0 | 1 | FMO_SHA256_00f691a1ddfb282c25ecce5c1a24902aeae3c7941ac54... |
| 112 | g2 | 0 | NaN | 0 | 0 | FMO_SHA256_0335e2f9610cbbc8189b4aa9c3698f91ee53fe6fe594f... |
What the additions changed¶
increments/<sample>/summary.json counts accepted peptides before (left, de novo only) and after (right, de novo plus molecular), and splits the gained peptides by whether they landed on an added candidate. protein_support.tsv.gz gives the same per protein.
In [5]:
Copied!
inc = json.loads((RUN / "increments/S01/summary.json").read_text())
keys = ["peptides_left", "peptides_right", "peptides_retained", "peptides_gained", "peptides_lost",
"left_candidates", "right_candidates", "candidate_sequences_added",
"one_peptide_gained_on_added_candidates", "one_peptide_gained_on_retained_candidates"]
pd.Series({k: inc[k] for k in keys}, name="S01")
inc = json.loads((RUN / "increments/S01/summary.json").read_text())
keys = ["peptides_left", "peptides_right", "peptides_retained", "peptides_gained", "peptides_lost",
"left_candidates", "right_candidates", "candidate_sequences_added",
"one_peptide_gained_on_added_candidates", "one_peptide_gained_on_retained_candidates"]
pd.Series({k: inc[k] for k in keys}, name="S01")
Out[5]:
peptides_left 548 peptides_right 550 peptides_retained 548 peptides_gained 2 peptides_lost 0 left_candidates 348 right_candidates 351 candidate_sequences_added 3 one_peptide_gained_on_added_candidates 2 one_peptide_gained_on_retained_candidates 0 Name: S01, dtype: int64
In [6]:
Copied!
support = pd.read_csv(RUN / "increments/S01/protein_support.tsv.gz", sep="\t")
support[support["gained_peptides"] > 0][["protein_sha256", "sequence_length", "in_left_candidates", "in_right_candidates", "left_peptides", "right_peptides", "gained_peptides"]]
support = pd.read_csv(RUN / "increments/S01/protein_support.tsv.gz", sep="\t")
support[support["gained_peptides"] > 0][["protein_sha256", "sequence_length", "in_left_candidates", "in_right_candidates", "left_peptides", "right_peptides", "gained_peptides"]]
Out[6]:
| protein_sha256 | sequence_length | in_left_candidates | in_right_candidates | left_peptides | right_peptides | gained_peptides | |
|---|---|---|---|---|---|---|---|
| 0 | 001f1586f7defaa3805fcb4cf73f6597963087434516ee2a929e8271... | 311 | 0 | 1 | 0 | 1 | 1 |
| 4 | 0335e2f9610cbbc8189b4aa9c3698f91ee53fe6fe594f78b52df78a2... | 137 | 0 | 1 | 0 | 1 | 1 |
Comparing the two searches¶
comparisons/<sample>/ holds the paired comparison of the two searches on the same spectra: the protein universe split into both, one side only and neither, the peptide overlap, and a per-ion LFQ comparison. When additions are few, the quantities of shared ions should not move; the log2 ratio distribution shows that.
In [7]:
Copied!
cmp = json.loads((RUN / "comparisons/S01/COMPLETE.json").read_text())
print("protein cells:", cmp["protein_cells"], "| universe:", cmp["protein_universe"])
lfq = pd.read_csv(RUN / "comparisons/S01/lfq_overlap.tsv", sep="\t")
paired = lfq.dropna(subset=["denovo_lfq", "denovo_plus_molecular_lfq"])
print(len(paired), "ions quantified in both searches; |log2 ratio| > 0.1 for", int((paired["log2_right_over_left"].abs() > 0.1).sum()))
fig, ax = plt.subplots(figsize=(5, 2.6))
ax.hist(paired["log2_right_over_left"], bins=41, color=OKABE_ITO[0])
ax.set_xlabel("log2 (de novo + molecular) / de novo, per ion"); ax.set_ylabel("ions")
plt.tight_layout()
cmp = json.loads((RUN / "comparisons/S01/COMPLETE.json").read_text())
print("protein cells:", cmp["protein_cells"], "| universe:", cmp["protein_universe"])
lfq = pd.read_csv(RUN / "comparisons/S01/lfq_overlap.tsv", sep="\t")
paired = lfq.dropna(subset=["denovo_lfq", "denovo_plus_molecular_lfq"])
print(len(paired), "ions quantified in both searches; |log2 ratio| > 0.1 for", int((paired["log2_right_over_left"].abs() > 0.1).sum()))
fig, ax = plt.subplots(figsize=(5, 2.6))
ax.hist(paired["log2_right_over_left"], bins=41, color=OKABE_ITO[0])
ax.set_xlabel("log2 (de novo + molecular) / de novo, per ion"); ax.set_ylabel("ions")
plt.tight_layout()
protein cells: {'both': 316, 'denovo_only': 0, 'denovo_plus_molecular_only': 2, 'neither': 33} | universe: 351
545 ions quantified in both searches; |log2 ratio| > 0.1 for 1
The same ledger and comparison files are produced for a real matched metagenome; the difference is only that the TPM values then come from read mapping. Molecular inputs describes the contract.