diff --git a/.gitignore b/.gitignore index b4b4e59..4356667 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,9 @@ evn* tmp .idea .coverage +*.egg-info +dist +build +test.ipynb +.tox +.python-version diff --git a/CHANGELOG.md b/CHANGELOG.md index 9bb100f..5a34708 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,23 @@ Fixes: - --> -## 0.0.13 (2025-03-31) +## 0.0.15 (2026-07-16) + +Features: + +- Replaced `ete3` runtime dependency with a custom `TreeNode`/`Tree` implementation in `src/pymutspec/annotation/phylo_tree.py` +- Added BioPython-based Newick parsing and an ete3-compatible tree API (node iteration, branch traversal, distance computation, node search, and Newick serialisation) +- Added `tests/test_tree_vs_ete3.py` to verify parity with ete3 `PhyloTree` for tree loading and node/edge iteration + +Fixes: + +- Kept `ete3` only in `dev` extras for comparison tests +- Fixed `get_tree_len` name to `get_tree_height` to better reflect its purpose (computing tree height from root to leaves) +- Removed requirement for trees have outgroup and be named "ROOT"; now any tree can be used, and the ingroup is defined as the sister node of the outgroup if present, or the root if not. + +**Full Changelog**: https://github.com/mitoclub/PyMutSpec/compare/0.0.14...0.0.15 + +## 0.0.14 (2025-03-31) Features: diff --git a/LICENSE b/LICENSE index b6f22dd..9f1ba85 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2025 mitofungen.com +Copyright (c) 2025 kpotoh Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index 540b720..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1 +0,0 @@ -include requirements.txt \ No newline at end of file diff --git a/README.md b/README.md index 48ffbc5..49f97af 100644 --- a/README.md +++ b/README.md @@ -9,65 +9,505 @@ https://pypi.org/project/pymutspec/) [![NeMu Paper](https://img.shields.io/badge/DOI-10.1093%2Fnar%2Fgkae438-blue)](https://doi.org/10.1093/nar/gkae438) - - +Python library for computing and visualising **mutational spectra** — genome-wide profiles of +single-nucleotide substitution patterns expressed in 12- or 192-component format (all possible +base changes with or without trinucleotide context). -Python library for mutational spectra analysis +PyMutSpec is the analysis layer used by the +[NeMu pipeline](https://nemu-pipeline.com) and has been applied to reconstruct neutral mutation +spectra from 2,591 chordate mitochondrial genomes +([Iliushchenko et al., bioRxiv 2023](https://doi.org/10.1101/2023.12.08.570826)). - +--- + +## Table of Contents + +1. [Installation](#installation) +2. [Development testing with tox](#development-testing-with-tox) +3. [Quick start](#quick-start) +4. [User guide](#user-guide) + - [Genetic codes and CodonAnnotation](#genetic-codes-and-codonannotation) + - [Computing expected mutation frequencies](#computing-expected-mutation-frequencies) + - [Observed mutations table format](#observed-mutations-table-format) + - [Calculating a mutational spectrum](#calculating-a-mutational-spectrum) + - [Plotting spectra](#plotting-spectra) + - [Collapsing a 192-component spectrum](#collapsing-a-192-component-spectrum) + - [Jackknife confidence intervals](#jackknife-confidence-intervals) + - [Edge-wise (per-branch) spectra](#edge-wise-per-branch-spectra) +5. [API overview](#api-overview) +6. [Links](#links) +7. [How to cite](#how-to-cite) + +--- ## Installation ```bash -pip3 install pymutspec +pip install pymutspec +``` + +## Development testing with tox + +Use tox to test the package in isolated environments across multiple Python versions. + +```bash +pyenv install 3.8 3.9 3.10 3.11 3.12 3.13 3.14 +pyenv local 3.8 3.9 3.10 3.11 3.12 3.13 3.14 + +pip install tox +tox p ``` -## Example code +Each tox environment installs the package and runs the test suite with pytest. + +If a specific Python interpreter is not installed locally, that environment will be skipped or fail to create. +Install missing versions with pyenv (or your system package manager) and re-run tox. + +--- + +## Quick start ```python +import pandas as pd from Bio import SeqIO -from pymutspec.annotation import calculate_mutspec, CodonAnnotation +from pymutspec.annotation import CodonAnnotation, calculate_mutspec from pymutspec.draw import plot_mutspec12, plot_mutspec192 -coda = CodonAnnotation(gencode=2) # mitochondrial genetic code +# 1. Initialise the codon annotator with the desired genetic code +# (2 = vertebrate mitochondrial; see NCBI for other codes) +coda = CodonAnnotation(gencode=2) + +# 2. Load the reference sequence for a single gene +ref_seq = str(next(SeqIO.parse("reference_gene.fasta", "fasta")).seq) + +# 3. Compute expected substitution frequencies from the reference +# Returns dicts keyed by label: 'all', 'syn', 'ff' (fourfold-degenerate) +sbs12_freqs, sbs192_freqs = coda.collect_exp_mut_freqs(ref_seq, labels=["syn"]) +exp_syn_12 = sbs12_freqs["syn"] # {sbs12: count, ...} e.g. {"C>A": 215.4, ...} +exp_syn_192 = sbs192_freqs["syn"] # {sbs192: count, ...} e.g. {"A[C>A]A": 8.2, ...} + +# 4. Load observed mutations (from NeMu-pipeline output or your own table) +obs = pd.read_csv("observed_mutations.tsv", sep="\t") + +# Filter to synonymous mutations only (Label == 1 or 2) +obs_syn = obs[obs["Label"] >= 1] + +# 5. Calculate spectra +ms12 = calculate_mutspec(obs_syn, exp_syn_12, use_context=False) +ms192 = calculate_mutspec(obs_syn, exp_syn_192, use_context=True) + +# 6. Plot +plot_mutspec12(ms12, title="12-component spectrum") +plot_mutspec192(ms192, title="192-component spectrum") +``` + +### Example output + + + + + +--- + +## User guide + +### Genetic codes and CodonAnnotation + +`CodonAnnotation` is the central class for working with codon-level mutation annotation. +It wraps a NCBI genetic code table and pre-computes synonymous and fourfold-degenerate codon +sets for fast lookup. + +```python +from pymutspec.annotation import CodonAnnotation + +# Standard genetic code (1) +coda_standard = CodonAnnotation(gencode=1) + +# Vertebrate mitochondrial code (2) — used for mtDNA analyses +coda_mito_vert = CodonAnnotation(gencode=2) + +# Inspect available methods +coda = CodonAnnotation(gencode=2) + +# Translate a codon +coda.translate_codon("ATG") # → 'M' + +# Check whether a mutation is synonymous +coda.is_syn_mut("CTA", "CTG") # → True (Leu → Leu) +coda.is_syn_mut("ATG", "ACG") # → False (Met → Thr) -path_to_observed_mutations = ... -path_to_reference_seq = ... +# Check whether a codon is fourfold-degenerate at the 3rd position +coda.is_fourfold("GTC") # → True (Val, all third-position changes are synonymous) +coda.is_fourfold("ATG") # → False -# load data (mutations and sequence) -gene = SeqIO.parse(path_to_reference_seq, format='fasta') -observed_mutations = pd.read_csv(path_to_observed_mutations, sep='\t') -for col in ['Mut', 'MutType']: - assert col in observed_mutations.columns +# Get mutation type label for a given codon pair and position-in-codon +label, aa_ref, aa_alt = coda.get_mut_type("CTA", "CTG", pic=2) +# label: 1 (syn), 2 (syn fourfold), 0 (non-syn), -1 (stop gain), ... +``` + +Available [NCBI genetic codes](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi): + +| Code | Description | +|------|-------------| +| 1 | Standard | +| 2 | Vertebrate mitochondrial | +| 4 | Mold, Protozoan, and Coelenterate mt | +| 5 | Invertebrate mitochondrial | +| 9 | Echinoderm and Flatworm mt | +| 13 | Ascidian mitochondrial | -# sample only syn mutations -mut_syn = observed_mutations[observed_mutations.MutType >= 1] # 0 for all mutations, 1 for syn, 2 for fourfold syn (syn4f) +--- -# derive expected mutations from reference gene -sbs12_freqs, sbs192_freqs = coda.collect_exp_mut_freqs(gene, labels['all', 'syn', 'syn4f']) -sbs12_freqs_syn = sbs12_freqs['syn'] -sbs192_freqs_syn = sbs192_freqs['syn'] +### Computing expected mutation frequencies -# calculate mutation spectra -spectra12 = calculate_mutspec(mut_syn, sbs12_freqs_syn, use_context=False) -spectra192 = calculate_mutspec(mut_syn, sbs192_freqs_syn, use_context=True) +Expected mutation frequencies capture how many opportunities each mutation type has in a given +sequence, accounting for codon structure and the requested synonymy class. + +```python +from pymutspec.annotation import CodonAnnotation -# plot mutation spectra -plot_mutspec12(spectra12) -plot_mutspec192(spectra192) +coda = CodonAnnotation(gencode=2) + +# From a single nucleotide sequence string (coding sequence, multiple of 3) +seq = "ATGCTAGTCGCATGA" # example coding sequence +sbs12_freqs, sbs192_freqs = coda.collect_exp_mut_freqs(seq, labels=["all", "syn", "ff"]) + +# sbs12_freqs["syn"] → dict mapping 12-component SBS to expected counts +# sbs192_freqs["syn"] → dict mapping 192-component SBS to expected counts + +# For a whole-genome / multi-gene analysis (e.g. from MSA or NeMu output), +# compute expected frequencies for each sequence in the alignment and average: +from Bio import SeqIO +sequences = [str(r.seq) for r in SeqIO.parse("msa_nuc.fasta", "fasta") + if r.id != "OUTGRP"] + +exp12_all, exp192_all = [], [] +for seq in sequences: + e12, e192 = coda.collect_exp_mut_freqs(seq, labels=["syn"]) + exp12_all.append(e12["syn"]) + exp192_all.append(e192["syn"]) + +# Average expected counts across sequences +import pandas as pd +exp12_mean = pd.DataFrame(exp12_all).mean().to_dict() +exp192_mean = pd.DataFrame(exp192_all).mean().to_dict() ``` -### Example spectra barplots +The `labels` parameter controls which mutation classes are included: + +| Label | Meaning | +|-------|---------| +| `"all"` | All single-nucleotide substitutions | +| `"syn"` | Synonymous substitutions only | +| `"ff"` | Fourfold-degenerate (most neutral) synonymous sites | + +--- - +### Observed mutations table format - +PyMutSpec works with observed-mutations tables produced by the +[NeMu pipeline](https://nemu-pipeline.com) or tables you construct yourself. +The minimum required columns are: + +| Column | Type | Description | +|--------|------|-------------| +| `Mut` | str | Mutation string in 192-component format: `X[N>M]Y` where `X`,`Y` are flanking nucleotides and `N>M` is the substitution, e.g. `A[C>T]G` | +| `Label` | int | Mutation type: `0` non-syn, `1` syn, `2` fourfold syn | + +Additional columns produced by NeMu and used by advanced functions: + +| Column | Description | +|--------|-------------| +| `ProbaFull` | Probability weight of the mutation (from ancestral reconstruction) | +| `RefNode` / `AltNode` | Parent and child node names in the phylogeny | +| `Site` | Position in gene (1-based) | +| `PosInCodon` | Position in codon (1-based: 1, 2, or 3) | +| `RefCodon` / `AltCodon` | Reference and mutated codons | +| `RefAa` / `AltAa` | Reference and mutated amino acids | + +```python +import pandas as pd + +obs = pd.read_csv("observed_mutations.tsv", sep="\t") +print(obs[["Mut", "Label", "ProbaFull"]].head()) +# Mut Label ProbaFull +# 0 A[C>T]G 1 0.91 +# 1 T[G>A]C 0 0.85 +# 2 C[A>G]T 2 0.78 +``` + +--- + +### Calculating a mutational spectrum + +`calculate_mutspec` divides observed counts by expected frequencies to produce a normalised +spectrum vector. + +```python +from pymutspec.annotation import calculate_mutspec + +# --- 12-component spectrum --- +ms12 = calculate_mutspec( + obs_muts=obs_syn, # DataFrame with at least "Mut" column + exp_muts=exp12_mean, # dict {sbs12: expected_count} + use_context=False, # False → 12-component + use_proba=False, # True → weight mutations by ProbaFull column + scale=True, # Normalise to sum to 1 +) +# ms12 columns: Mut, ObsNum, ExpNum, MutSpec + +# --- 192-component spectrum --- +ms192 = calculate_mutspec( + obs_muts=obs_syn, + exp_muts=exp192_mean, # dict {sbs192: expected_count} + use_context=True, # True → 192-component + use_proba=True, # Use ancestral-reconstruction probabilities + scale=True, +) +``` + +When using NeMu-pipeline output with probability weights, set `use_proba=True` and make sure +the `ProbaFull` column is present. + +The `Label` column codes for mutation types: + +| Value | Meaning | +|-------|---------| +| `0` | Non-synonymous | +| `1` | Synonymous | +| `2` | Fourfold-degenerate synonymous | +| `-1` | Stop gain | +| `-2` | Stop loss | +| `-3` | Stop-to-stop | + +Filter to the mutation class you need before calling `calculate_mutspec`: + +```python +# Synonymous only (Label == 1 or 2) +obs_syn = obs[obs["Label"] >= 1] + +# Fourfold-degenerate synonymous only +obs_ff = obs[obs["Label"] == 2] + +# All mutations +obs_all = obs.copy() +``` + +--- + +### Plotting spectra + +```python +from pymutspec.draw import plot_mutspec12, plot_mutspec192 + +# 12-component barplot +plot_mutspec12(ms12, title="Human CytB – synonymous spectrum", ylabel="MutSpec") + +# 192-component barplot (COSMIC order, default) +plot_mutspec192(ms192, title="192-component spectrum") + +# 192-component barplot with KK-style compact labels +from pymutspec.draw.spectra import ordered_sbs192_kk +plot_mutspec192(ms192, sbs_order=ordered_sbs192_kk, labels_style="kk", + title="192-component spectrum (KK order)") + +# Save to file without displaying +plot_mutspec192(ms192, savepath="spectrum.png", show=False) + +# Draw on an existing axes (useful for multi-panel figures) +import matplotlib.pyplot as plt +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4)) +plot_mutspec12(ms12, ax=ax1, show=False, title="12-comp") +plot_mutspec12(ms12_2, ax=ax2, show=False, title="Other sample") +plt.tight_layout() +plt.show() +``` + +`labels_style` options for 192-component plots: + +| Value | Tick label format | Example | +|------------|-------------------|---------| +| `"cosmic"` | Full COSMIC string | `A[C>A]T` | +| `"long"` | Same as `"cosmic"` | `A[C>A]T` | +| `"kk"` | Compact KK format | `CA: ACT` | + +--- + +### Collapsing a 192-component spectrum + +When you want to aggregate per-branch or per-sample 192-component spectra into a single matrix +and then reduce to 12 components: + +```python +from pymutspec.annotation import complete_sbs192_columns +from pymutspec.annotation.spectra import collapse_sbs192 +from pymutspec.constants import possible_sbs192 + +# Suppose you have a list of per-sample spectra as dicts +samples = [ + {"A[C>A]A": 0.012, "T[G>T]T": 0.008, ...}, # sample 1 + {"A[C>A]A": 0.009, "C[C>T]G": 0.021, ...}, # sample 2 +] +df192 = pd.DataFrame(samples) + +# Ensure all 192 columns are present (fill missing with 0) +df192 = complete_sbs192_columns(df192) # shape: (n_samples, 192) + +# Collapse to 12 components +df12 = collapse_sbs192(df192, to=12) # shape: (n_samples, 12) +``` + +--- + +### Jackknife confidence intervals + +`jackknife_spectra_sampling` estimates spectrum variability by repeatedly sampling a fraction +of branches (tree edges) without replacement: + +```python +from pymutspec.annotation.spectra import jackknife_spectra_sampling + +# obs: observed mutations DataFrame (long format with RefNode, AltNode, Mut, ProbaFull) +# exp: expected frequencies DataFrame (long format with Node, Mut, Proba) +spectra_samples = jackknife_spectra_sampling(obs, exp, frac=0.5, n=1000) +# spectra_samples: DataFrame of shape (1000, 192) + +# Derive confidence intervals +q05 = spectra_samples.quantile(0.05) +q95 = spectra_samples.quantile(0.95) +median = spectra_samples.median() + +# Build a summary table for plotting with error bars +import pandas as pd +summary = pd.DataFrame({ + "Mut": possible_sbs192, + "MutSpec": ms192.set_index("Mut")["MutSpec"], + "MutSpec_median": median.values, + "MutSpec_q05": q05.values, + "MutSpec_q95": q95.values, +}) +# plot_mutspec192 will automatically draw error bars if those columns are present +plot_mutspec192(summary, title="Spectrum with 90% CI") +``` + +--- + +### Edge-wise (per-branch) spectra + +For phylogenomic analyses it is useful to compute a separate mutational spectrum for each +branch of the phylogeny, then compare them: + +```python +from pymutspec.annotation.spectra import calc_edgewise_spectra, get_cossim + +spectra_per_edge = calc_edgewise_spectra( + obs, # observed mutations with RefNode/AltNode columns + exp, # expected frequencies with Node column + nmtypes_cutoff=10, # minimum distinct mutation types per branch + nobs_cuttof=10, # minimum observed mutations per branch + scale=True, # normalise each branch spectrum +) +# spectra_per_edge: DataFrame indexed by (RefNode, AltNode), 192 columns + +# Compare two sets of branch spectra using cosine similarity +cossim = get_cossim(spectra_per_edge.loc[["nodeA"]], spectra_per_edge.loc[["nodeB"]]) +``` + +--- + +## API overview + +### `pymutspec.annotation` + +| Function / Class | Description | +|------------------|-------------| +| `CodonAnnotation(gencode)` | Codon-level annotator for a given NCBI genetic code | +| `CodonAnnotation.collect_exp_mut_freqs(seq, labels)` | Expected substitution counts from a nucleotide sequence | +| `CodonAnnotation.collect_exp_muts(seq, labels)` | Expected substitutions as a long-format DataFrame | +| `CodonAnnotation.is_syn_mut(cdn1, cdn2)` | Test if a codon change is synonymous | +| `CodonAnnotation.is_fourfold(cdn)` | Test if a codon is fourfold-degenerate | +| `CodonAnnotation.translate_codon(cdn)` | Translate a codon to a single-letter amino acid | +| `CodonAnnotation.get_mut_type(cdn1, cdn2, pic)` | Return mutation type label, ref AA, alt AA | +| `calculate_mutspec(obs, exp, ...)` | Compute a 12- or 192-component mutational spectrum | +| `complete_sbs192_columns(df)` | Ensure a DataFrame has all 192 SBS columns (fill missing with 0) | +| `mutations_summary(mutations, ...)` | Tabulate mutation type counts from an observed-mutations table | +| `rev_comp(sbs)` | Reverse-complement a 192-component SBS string, e.g. `A[C>A]T` → `A[G>T]T` | +| `lbl2lbl_id(lbl)` / `lbl_id2lbl(id)` | Convert between label string and integer code | + +### `pymutspec.annotation.spectra` + +| Function | Description | +|----------|-------------| +| `collapse_sbs192(df, to=12)` | Sum 192-component spectra to 12 components | +| `collapse_mutspec(ms192)` | Collapse a 192-row spectrum DataFrame to 96 components via strand symmetry | +| `filter_outlier_branches(obs_df, ...)` | Remove branches with outlier-high mutation counts (IQR method) | +| `jackknife_spectra_sampling(obs, exp, frac, n)` | Bootstrap spectrum confidence intervals by branch resampling | +| `calc_edgewise_spectra(obs, exp, ...)` | Per-branch (edge-wise) mutational spectra | +| `get_cossim(a, b)` | Row-wise cosine similarity between two spectrum DataFrames | +| `get_eucdist(a, b)` | Row-wise Euclidean distance between two spectrum DataFrames | + +### `pymutspec.draw` + +| Function | Description | +|----------|-------------| +| `plot_mutspec(mutspec, sbs_kind, ...)` | Generic plotting function (12- or 192-component) | +| `plot_mutspec12(mutspec, ...)` | Convenience wrapper for 12-component barplots | +| `plot_mutspec192(mutspec, ...)` | Convenience wrapper for 192-component barplots | + +--- ## Links -1. [IQ-Tree2](http://www.iqtree.org/) - efficient software for phylogenomic inference -2. [Genetic codes](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi?chapter=tgencodes#SG1) +1. [NeMu pipeline](https://nemu-pipeline.com) — end-to-end pipeline for neutral mutation spectra from evolutionary data +2. [IQ-Tree2](http://www.iqtree.org/) — efficient software for phylogenomic inference +3. [Genetic codes](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi) — NCBI codon tables +4. [mtDNA chordata analysis](https://github.com/mitoclub/mtdna-192component-mutspec-chordata) — example notebooks using PyMutSpec on 2,591 vertebrate mitochondrial genomes + +--- + +## How to cite + +If you use PyMutSpec in your work, please cite the paper that describes the methods: + +Efimenko, B., Popadin, K., & Gunbin, K. (2024). NeMu: a comprehensive pipeline for accurate +reconstruction of neutral mutation spectra from evolutionary data. *Nucleic Acids Research*, +52(W1), W108–W115. + +Suggested BibTeX entry: + +```bibtex +@article{Efimenko2024NeMu, + author = {Efimenko, Bogdan and Popadin, Konstantin and Gunbin, Konstantin}, + title = {NeMu: a comprehensive pipeline for accurate reconstruction of neutral mutation spectra from evolutionary data}, + journal = {Nucleic Acids Research}, + volume = {52}, + number = {W1}, + pages = {W108--W115}, + year = {2024}, + doi = {10.1093/nar/gkae438}, +} +``` + +Thank you for citing the work if PyMutSpec aids your research. + +--- + + -## How to cite? +## TODO - Bogdan Efimenko, Konstantin Popadin, Konstantin Gunbin, NeMu: a comprehensive pipeline for accurate reconstruction of neutral mutation spectra from evolutionary data, Nucleic Acids Research, Volume 52, Issue W1, 5 July 2024, Pages W108–W115, https://doi.org/10.1093/nar/gkae438 +- [x] Custom tree implementation +- [ ] New way of annotation from HGT +- [ ] Separate scripts for mutations collection and annotation +- [ ] Integrate parallelisation from HGT project +- [ ] Add feature to calculate mutation rate (MutRate) +- [ ] Rename some functions and variables for better readability diff --git a/pymutspec/__init__.py b/pymutspec/__init__.py deleted file mode 100644 index 31898f5..0000000 --- a/pymutspec/__init__.py +++ /dev/null @@ -1 +0,0 @@ -__version__ = "0.0.13" \ No newline at end of file diff --git a/pymutspec/annotation/auxiliary.py b/pymutspec/annotation/auxiliary.py deleted file mode 100644 index 5df2ae8..0000000 --- a/pymutspec/annotation/auxiliary.py +++ /dev/null @@ -1,31 +0,0 @@ -transcriptor = str.maketrans("ACGT", "TGCA") - - -def rev_comp(mut: str): - new_mut = mut[-1] + mut[1:-1] + mut[0] - new_mut = new_mut.translate(transcriptor) - return new_mut - - -def lbl_id2lbl(lbl_id: int) -> str: - if lbl_id == 0: - lbl = "all" - elif lbl_id == 1: - lbl = "syn" - elif lbl_id == 2: - lbl = "ff" - else: - raise NotImplementedError() - return lbl - - -def lbl2lbl_id(lbl: str) -> int: - if lbl == "all": - lbl_id = 0 - elif lbl == "syn" or lbl == "syn_c": - lbl_id = 1 - elif lbl == "ff" or lbl == "syn4f": - lbl_id = 2 - else: - raise NotImplementedError() - return lbl_id diff --git a/pymutspec/annotation/tree.py b/pymutspec/annotation/tree.py deleted file mode 100644 index 556958b..0000000 --- a/pymutspec/annotation/tree.py +++ /dev/null @@ -1,82 +0,0 @@ -from queue import Queue -from statistics import geometric_mean - -import numpy as np -from ete3 import PhyloTree, PhyloNode - - -def node_parent(node: PhyloNode): - try: - return next(node.iter_ancestors()) - except BaseException: - return None - - -def iter_tree_edges(tree: PhyloTree): - discovered_nodes = set() - discovered_nodes.add(tree.name) - Q = Queue() - Q.put(tree) - - while not Q.empty(): - cur_node = Q.get() - for child in cur_node.children: - Q.put(child) - - if cur_node.name not in discovered_nodes: - discovered_nodes.add(cur_node.name) - alt_node = cur_node - ref_node = node_parent(alt_node) - yield ref_node, alt_node - - -def get_tree_len(tree: PhyloTree, mode='geom_mean'): - ''' - TODO check if tree is rooted - - Params: - - mode: str - calculate 'mean', 'geom_mean' or 'max' of distribution of len from current node to leaves - ''' - assert tree.name != 'ROOT' - - if mode == 'max': - _, md = tree.get_farthest_leaf() - elif mode in ['mean', 'geom_mean']: - distances_to_leaves = [] - for leaf in tree.iter_leaves(): - d = tree.get_distance(leaf) - distances_to_leaves.append(d) - - if mode == 'mean': - md = np.mean(distances_to_leaves) - elif mode == 'geom_mean': - md = geometric_mean(distances_to_leaves) - - else: - raise TypeError(f"mode must be 'mean', 'geom_mean' or 'max'") - - return md - - -def get_ingroup_root(tree: PhyloTree) -> PhyloTree: - assert len(tree.children) == 2, 'Tree must be binary' - found_outgroup = False - for node in tree.children: - if node.is_leaf(): - found_outgroup = True - else: - ingrp = node - - if found_outgroup: - return ingrp - else: - return tree - - -def calc_phylocoefs(tree: PhyloTree): - tree_len = get_tree_len(get_ingroup_root(tree), 'geom_mean') - phylocoefs = {tree.name: 1 - min(0.999, tree.get_closest_leaf()[1] / tree_len)} - for node in tree.iter_descendants(): - _closest, d = node.get_closest_leaf() - phylocoefs[node.name] = 1 - min(0.99999, d / tree_len) - return phylocoefs diff --git a/pymutspec/draw/__init__.py b/pymutspec/draw/__init__.py deleted file mode 100644 index e980a0f..0000000 --- a/pymutspec/draw/__init__.py +++ /dev/null @@ -1,4 +0,0 @@ -from .spectra import plot_mutspec12, plot_mutspec192, plot_mutspec192kk, _prepare_nice_labels, plot_mutspec96 - -def plot_mutspec192box(*args, **kwargs): - print("WARNING: this alias function is removed! Use plot_mutspec192(style='box') instead") diff --git a/pymutspec/draw/sbs_orders.py b/pymutspec/draw/sbs_orders.py deleted file mode 100644 index c5ed8bc..0000000 --- a/pymutspec/draw/sbs_orders.py +++ /dev/null @@ -1,22 +0,0 @@ -import pandas as pd - -from ..constants import possible_sbs192 -from ..annotation import rev_comp, transcriptor - -kk_lbls = "A>C A>G A>T C>T G>C G>T".split() -cosmic_lbls = "C>A C>G C>T T>A T>C T>G".split() - -df = pd.DataFrame({"sbs": possible_sbs192}) -df["sbs_base"] = df["sbs"].str.slice(2, 5) -df["sbs_base_revcomp"] = df["sbs_base"].str.translate(transcriptor) -df["sbs_revcomp"] = df["sbs"].apply(rev_comp) -df["is_cosmic"] = df["sbs_base"].isin(cosmic_lbls) -df["is_kk"] = df["sbs_base"].isin(kk_lbls) -df["sbs_base_for_sorting_kp"] = df.apply( - lambda x: x.sbs_base + "1" if x.is_cosmic else x.sbs_base_revcomp + "2", axis=1) -df["sbs_for_ordering_kk"] = df.apply(lambda x: x.sbs if x.is_kk else x.sbs_revcomp, axis=1) -df["sbs_for_ordering_kp"] = df.apply(lambda x: x.sbs if x.is_cosmic else x.sbs_revcomp, axis=1) - -ordered_sbs192_kp = list(df.sort_values(["sbs_base_for_sorting_kp", "sbs_for_ordering_kp"]).sbs.values) -ordered_sbs192_kk = list(df.sort_values(["sbs_base", "sbs_for_ordering_kk"]).sbs.values) -del df diff --git a/pymutspec/draw/spectra.py b/pymutspec/draw/spectra.py deleted file mode 100644 index cd4bc43..0000000 --- a/pymutspec/draw/spectra.py +++ /dev/null @@ -1,362 +0,0 @@ -""" -Functionality to plot mutational spectrums -""" - -from typing import Iterable - -import pandas as pd -import matplotlib -import matplotlib.pyplot as plt -import seaborn as sns - -from ..constants import possible_sbs192, possible_sbs96 -from .sbs_orders import ordered_sbs192_kk, ordered_sbs192_kp - -color_mapping6 = { - "C>A": "deepskyblue", - "C>G": "black", - "C>T": "red", - "T>A": "silver", - "T>C": "green", - "T>G": "pink", -} -color_mapping12 = { - "C>A": "deepskyblue", - "G>T": "deepskyblue", - "C>G": "black", - "G>C": "black", - "C>T": "red", - "G>A": "red", - "T>A": "silver", - "A>T": "silver", - "T>C": "yellowgreen", - "A>G": "yellowgreen", - "T>G": "pink", - "A>C": "pink", -} -sbs12_ordered = ["C>A", "G>T", "C>G", "G>C", "C>T", "G>A", "T>A", "A>T", "T>C", "A>G", "T>G", "A>C"] -_smpl = pd.DataFrame({"Mut": possible_sbs192}) -_smpl["MutBase"] = _smpl.Mut.str.slice(2, 5) -_smpl["Context"] = _smpl.Mut.str.get(0) + _smpl.Mut.str.get(2) + _smpl.Mut.str.get(-1) -colors192 = _smpl.sort_values(["MutBase", "Context"])["MutBase"].map(color_mapping12).values -colors12 = [color_mapping12[sbs] for sbs in sbs12_ordered] - - -def _prepare_nice_labels(sbs192: Iterable[str], kk=False): - _nice_sbs = [] - prev = None - for sbs in sbs192: - if prev is not None and sbs[2:5] != prev[2:5]: - _nice_sbs.append("") - sbs_nice = sbs[2] + sbs[4] + ": " + sbs[0] + sbs[2] + sbs[-1] if kk else sbs - _nice_sbs.append(sbs_nice) - prev = sbs - return _nice_sbs - - -def _coloring192kk(): - colors = "red yellow lime blue".split() - while True: - for clr in colors: - yield clr - - -def plot_mutspec12( - mutspec: pd.DataFrame, - spectra_col="MutSpec", - title="Full mutational spectrum", - ylabel=None, - figsize=(6, 4), - style="bar", - savepath=None, - fontname=None, - ticksize=8, - titlesize=14, - ylabelsize=12, - show=True, - ax=None, - **kwargs, - ): - # TODO add checks of mutspec12 - # TODO add description to all plot* functions - - if ax is None: - fig = plt.figure(figsize=figsize) - ax = fig.gca() - elif isinstance(ax, matplotlib.axes._subplots.AxesSubplot): - pass - else: - raise ValueError('ax must be None or matplotlib.axes._subplots.AxesSubplot') - - if style == "bar": - _cols = set(mutspec.columns) - if 'MutSpec_median' in _cols and 'MutSpec_q05' in _cols and 'MutSpec_q95' in _cols: - ax = sns.barplot(data=mutspec, x="Mut", y='MutSpec', - order=sbs12_ordered, ax=ax, **kwargs) - mutspec_ordered = mutspec.set_index('Mut').loc[sbs12_ordered] - mutspec_ordered['MutSpec_q05'] = mutspec_ordered['MutSpec_median'] - mutspec_ordered['MutSpec_q05'] - mutspec_ordered['MutSpec_q95'] -= mutspec_ordered['MutSpec_median'] - ax.errorbar(mutspec_ordered.index, mutspec_ordered['MutSpec_median'], - yerr=mutspec_ordered[['MutSpec_q05', 'MutSpec_q95']].values.T, - fmt=".", color="gray", elinewidth=0.7, capsize=4) - else: - ax = sns.barplot(data=mutspec, x="Mut", y=spectra_col, - order=sbs12_ordered, ax=ax, **kwargs) - - elif style == "box": - ax = sns.boxplot(x="Mut", y=spectra_col, data=mutspec, - order=sbs12_ordered, ax=ax, **kwargs) - - else: - raise NotImplementedError - - ax.grid(axis="y", alpha=.7, linewidth=0.5) - - # map colors to bars - for bar, clr in zip(ax.patches, colors12): - bar.set_color(clr) - - ax.set_title(title, fontsize=titlesize, fontname=fontname) - ax.set_ylabel(ylabel if ylabel else "", fontsize=ylabelsize, fontname=fontname) - ax.set_xlabel("") - - plt.xticks(fontsize=ticksize, fontname=fontname) - - if savepath is not None: - plt.savefig(savepath, bbox_inches="tight") - if show: - plt.show() - else: - plt.close() - return ax - - -def plot_mutspec192( - mutspec192: pd.DataFrame, - spectra_col="MutSpec", - title="Mutational spectrum", - ylabel=None, - figsize=(24, 8), - style="bar", - labels_style="cosmic", - sbs_order=ordered_sbs192_kp, - savepath=None, - fontname=None, - ticksize=6, - titlesize=16, - ylabelsize=16, - show=True, - ax=None, - **kwargs, - ): - """ - Plot barblot of given mutational spectrum calculated from single nucleotide substitutions - - Arguments - --------- - mutspec192: pd.DataFrame - table, containing 192 component mutational spectrum for one or many species, all substitutions must be presented in the table - title: str, default = 'Mutational spectrum' - Title on the plot - savepath: str, default = None - Path to output plot file. If None no images will be written - labels_style: str, default = 'cosmic' - 'cosmic': A[C>T]T, 'long': CT: ACT - """ - if "filepath" in kwargs: - savepath = kwargs["filepath"] - print("savepath =", savepath) - kwargs.pop("filepath") - - # TODO add checks of mutspec192 - ms192 = mutspec192.copy() - if labels_style == "long": - ms192["long_lbl"] = ms192.Mut.str.get(2) + ms192.Mut.str.get(4) + ": " + ms192.Mut.str.get(0) + ms192.Mut.str.get(2) + ms192.Mut.str.get(-1) - order = _prepare_nice_labels(sbs_order, kk=True) - x_col = "long_lbl" - elif labels_style == "cosmic": - order = _prepare_nice_labels(sbs_order, kk=False) - x_col = "Mut" - else: - raise ValueError("Available labels_style are: 'cosmic' and 'long'") - - if ax is None: - fig = plt.figure(figsize=figsize) - ax = fig.gca() - elif isinstance(ax, matplotlib.axes._subplots.AxesSubplot): - pass - else: - raise ValueError('ax must be None or matplotlib.axes._subplots.AxesSubplot') - - if style == "bar": - _cols = set(ms192.columns) - if 'MutSpec_median' in _cols and 'MutSpec_q05' in _cols and 'MutSpec_q95' in _cols: - ax = sns.barplot(data=ms192, x='Mut', y='MutSpec', order=order, ax=ax, **kwargs) - mutspec_ordered = ms192.set_index('Mut').loc[sbs_order] - mutspec_ordered['MutSpec_q05'] = mutspec_ordered['MutSpec_median'] - mutspec_ordered['MutSpec_q05'] - mutspec_ordered['MutSpec_q95'] -= mutspec_ordered['MutSpec_median'] - - ymed, yerr_min, yerr_max = [], [], [] - for mt in order: - if mt == '': - ymed.append(0.) - yerr_min.append(0.) - yerr_max.append(0.) - else: - ymed.append(mutspec_ordered.loc[mt, 'MutSpec_median']) - yerr_min.append(mutspec_ordered.loc[mt, 'MutSpec_q05']) - yerr_max.append(mutspec_ordered.loc[mt, 'MutSpec_q95']) - yerr = [yerr_min, yerr_max] - x = list(range(len(order))) - ax.errorbar(x, ymed, yerr=yerr, fmt=".", color="gray", elinewidth=0.7, capsize=2) - else: - ax = sns.barplot(data=ms192, x=x_col, y=spectra_col, order=order, errwidth=1, ax=ax, **kwargs) - elif style == "box": - ax = sns.boxplot( - x=x_col, y=spectra_col, data=ms192, order=order, ax=fig.gca(), **kwargs, - ) - # plt.savefig(savepath, dpi=300, bbox_inches="tight") - # return - ax.grid(axis="y", alpha=.7, linewidth=0.5) - ax.set_title(title, fontsize=titlesize, fontname=fontname) - ax.set_ylabel(ylabel if ylabel else "", fontsize=ylabelsize, fontname=fontname) - ax.set_xlabel("") - # map colors to bars - width = 0.4 - shift = None - for bar, sbs in zip(ax.patches, order): - if len(sbs): - s = sbs[0] + ">" + sbs[1] if labels_style == "long" else sbs[2:5] - bar.set_color(color_mapping12[s]) - bar.set_alpha(alpha=0.9) - if style == "bar": - if not shift: - # calculate one time instead of 192 - shift = (bar.get_width() - width) / 2 - bar.set_width(width) - bar.set_x(bar.get_x() + shift) - - plt.xticks(rotation=90, fontsize=ticksize, fontname=fontname) - - if savepath is not None: - plt.savefig(savepath, dpi=300, bbox_inches="tight") - if show: - plt.show() - else: - plt.close() - return ax - - -def plot_mutspec96( - spectra: pd.DataFrame, - ylabel="MutSpec", - title="Mutational spectrum", - figsize=(15, 5), - style="bar", - labels_style="cosmic", - sbs_order=possible_sbs96, - savepath=None, - fontsize=9, - titlesize=16, - fontname=None, - show=True, - **kwargs, - ): - """ - Plot barblot of given mutational spectrum calculated from single nucleotide substitutions - - Arguments - --------- - spectra: pd.DataFrame - table, containing 192 component mutational spectrum for one or many species, all substitutions must be presented in the table - title: str, default = 'Mutational spectrum' - Title on the plot - savepath: str, default = None - Path to output plot file. If None no images will be written - labels_style: str, default = 'cosmic' - 'cosmic': A[C>T]T, 'long': CT: ACT - """ - if "filepath" in kwargs: - savepath = kwargs["filepath"] - print("savepath =", savepath) - kwargs.pop("filepath") - - # TODO add checks of spectra - ms96 = spectra.copy() - if labels_style == "long": - ms96["long_lbl"] = ms96.Mut.str.get(2) + ms96.Mut.str.get(4) + ": " + ms96.Mut.str.get(0) + ms96.Mut.str.get(2) + ms96.Mut.str.get(-1) - order = _prepare_nice_labels(sbs_order, kk=True) - x_col = "long_lbl" - elif labels_style == "cosmic": - order = _prepare_nice_labels(sbs_order, kk=False) - x_col = "Mut" - else: - raise ValueError("Available labels_style are: 'cosmic' and 'long'") - fig = plt.figure(figsize=figsize) - if style == "bar": - ax = sns.barplot( - x=x_col, y=ylabel, data=ms96, order=order, errwidth=1, ax=fig.gca(), **kwargs, - ) - elif style == "box": - ax = sns.boxplot( - x=x_col, y=ylabel, data=ms96, order=order, ax=fig.gca(), **kwargs, - ) - ax.grid(axis="y", alpha=.7, linewidth=0.5) - ax.set_title(title, fontsize=titlesize, fontname=fontname) - ax.set_xlabel("") - ax.set_ylabel("") - # map colors to bars - width = 0.4 - shift = None - for bar, sbs in zip(ax.patches, order): - if len(sbs): - s = sbs[0] + ">" + sbs[1] if labels_style == "long" else sbs[2:5] - bar.set_color(color_mapping6[s]) - bar.set_alpha(alpha=0.9) - if style == "bar": - if not shift: - # calculate one time instead of 192 - shift = (bar.get_width() - width) / 2 - bar.set_width(width) - bar.set_x(bar.get_x() + shift) - - plt.xticks(rotation=90, fontsize=fontsize, fontname=fontname) - - if savepath is not None: - plt.savefig(savepath, dpi=300, bbox_inches="tight") - if show: - plt.show() - else: - plt.close() - return ax - - -def plot_mutspec192kk(mutspec192: pd.DataFrame, ylabel="MutSpec", title="Mutational spectrum", show=True, figsize=(24, 6), filepath=None): - ms192 = mutspec192.copy() - ms192["long_lbl"] = ms192.Mut.str.get(2) + ms192.Mut.str.get(4) + ": " + ms192.Mut.str.get(0) + ms192.Mut.str.get(2) + ms192.Mut.str.get(-1) - fig = plt.figure(figsize=figsize) - ax = fig.add_subplot(111) - ax.grid(axis="y", alpha=.7, linewidth=0.5) - order = _prepare_nice_labels(ordered_sbs192_kk, True) - sns.barplot( - x="long_lbl", y=ylabel, data=ms192, - order=order, - errwidth=1, ax=fig.gca(), - ) - plt.xticks(rotation=90, fontsize=7, fontname=None) - ax.set_title(title) - ax.set_xlabel("") - ax.set_ylabel("Mutational spectrum") - # map colors to bars - clrs_iterator = _coloring192kk() - for bar, sbs in zip(ax.patches, order): - if len(sbs): - bar.set_color(next(clrs_iterator)) - bar.set_alpha(alpha=0.9) - bar.set_width(0.3) - if filepath is not None: - plt.savefig(filepath, dpi=300, bbox_inches="tight") - if show: - plt.show() - else: - plt.close() diff --git a/pyproject.toml b/pyproject.toml index 196cb89..b2d676b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,10 @@ classifiers = [ "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", "Topic :: Scientific/Engineering", "Topic :: Scientific/Engineering :: Bio-Informatics", "Topic :: Software Development :: Libraries :: Python Modules", @@ -26,27 +30,26 @@ classifiers = [ ] requires-python = ">=3.8" dependencies = [ - "numpy<2.0", - "pandas<2.0", - "seaborn==0.12.2", + "numpy", + "pandas", + "seaborn", "biopython>=1.75", - "ete3>=3", - "tqdm", "PyYAML", - "pytest", "click", ] [project.optional-dependencies] -dev = ["pytest", "pytest-cov", "flake8"] +dev = ["pytest", "pytest-cov", "flake8", "ete3>=3", "legacy-cgi", "tox", "ruff", "build", "twine"] [project.urls] Homepage = "https://github.com/mitoclub/PyMutSpec" Changelog = "https://github.com/mitoclub/PyMutSpec/blob/master/CHANGELOG.md" "Bug Tracker" = "https://github.com/mitoclub/PyMutSpec/issues" +[tool.setuptools.packages.find] +where = ["src"] + [tool.setuptools] -packages = ["pymutspec"] script-files = ["scripts/multifasta_coding.py", "scripts/alignment2iqtree_states.py", "scripts/select_records.py", "scripts/iqtree_states_add_part.py", "scripts/collect_mutations.py", "scripts/calculate_mutspec.py"] @@ -54,3 +57,7 @@ script-files = ["scripts/multifasta_coding.py", "scripts/alignment2iqtree_states [tool.setuptools.dynamic] version = {attr = "pymutspec.__version__"} # any module attribute compatible with ast.literal_eval + +[tool.pytest.ini_options] +pythonpath = ["src"] +testpaths = "tests" diff --git a/requirements.dev.txt b/requirements.dev.txt deleted file mode 100644 index eabcf45..0000000 --- a/requirements.dev.txt +++ /dev/null @@ -1,15 +0,0 @@ -numpy<2.0 -pandas<2.0 -matplotlib~=3.5.1 -seaborn==0.12.2 -biopython>=1.75 -ete3>=3 -tqdm -PyYAML -pytest -click -# scipy~=1.8.1 -# statsmodels~=0.13.2 -pytest -pytest-cov -flake8 \ No newline at end of file diff --git a/scripts/1.terminal_genomes2iqtree_format.py b/scripts/1.terminal_genomes2iqtree_format.py index 1e2b002..cef296b 100644 --- a/scripts/1.terminal_genomes2iqtree_format.py +++ b/scripts/1.terminal_genomes2iqtree_format.py @@ -29,7 +29,7 @@ def parse_alignment_and_write_states(files: list, outfile) -> Tuple[str, int]: return states table and full alignment length """ - print(f"Processing...", file=sys.stderr) + print("Processing...", file=sys.stderr) handle = open(outfile, "w") ngenes = len(files) columns = "Node Part Site State p_A p_C p_G p_T".split() diff --git a/scripts/alignment2iqtree_states.py b/scripts/alignment2iqtree_states.py index 9a48f03..7d63a98 100644 --- a/scripts/alignment2iqtree_states.py +++ b/scripts/alignment2iqtree_states.py @@ -24,7 +24,7 @@ def parse_alignment_and_write_states(files: list, outfile) -> Tuple[str, int]: return states table and full alignment length """ - print(f"Processing...", file=sys.stderr) + print("Processing...", file=sys.stderr) handle = open(outfile, "w") ngenes = len(files) columns = "Node Part Site State p_A p_C p_G p_T".split() @@ -72,7 +72,7 @@ def parse_alignment_and_write_states(files: list, outfile) -> Tuple[str, int]: @click.argument("states", nargs=1, type=click.Path(writable=True)) def main(alignment, states): parse_alignment_and_write_states(alignment, states) - print(f"Done.", file=sys.stderr) + print("Done.", file=sys.stderr) if __name__ == "__main__": main() diff --git a/scripts/calculate_mutspec.py b/scripts/calculate_mutspec.py index 3e606e8..3742087 100644 --- a/scripts/calculate_mutspec.py +++ b/scripts/calculate_mutspec.py @@ -2,7 +2,6 @@ import os import sys -from functools import partial import click import pandas as pd @@ -12,7 +11,6 @@ ) from pymutspec.constants import possible_sbs12, possible_sbs192 from pymutspec.draw import plot_mutspec12, plot_mutspec192 -from pymutspec.draw.sbs_orders import ordered_sbs192_kp def save_tsv(df: pd.DataFrame, path): diff --git a/scripts/collect_mutations.py b/scripts/collect_mutations.py index a8b9a27..0241dc6 100644 --- a/scripts/collect_mutations.py +++ b/scripts/collect_mutations.py @@ -13,7 +13,7 @@ import click import numpy as np import pandas as pd -from ete3 import PhyloTree +from pymutspec.annotation.phylo_tree import Tree from pymutspec.annotation import ( CodonAnnotation, calculate_mutspec, @@ -77,7 +77,7 @@ def __init__( self.mut_labels.append("nonsyn") logger.info(f"Types of mutations to collect and process: {self.mut_labels}") self.fp_format = np.float32 - self.tree = PhyloTree(path_to_tree, format=1) + self.tree = Tree(path_to_tree, format=1) logger.info( f"Tree loaded, number of leaf nodes: {len(self.tree)}, " f"total number of nodes: {len(self.tree.get_cached_content())}" diff --git a/scripts/collect_mutations_parallel.py b/scripts/collect_mutations_parallel.py index 2f202c3..30c3f3f 100644 --- a/scripts/collect_mutations_parallel.py +++ b/scripts/collect_mutations_parallel.py @@ -17,9 +17,9 @@ import numpy as np import pandas as pd -from ete3 import PhyloTree +from pymutspec.annotation.phylo_tree import Tree from pymutspec.annotation import ( - CodonAnnotation, iter_tree_edges, get_tree_len, + CodonAnnotation, iter_tree_edges, get_tree_height, ) from pymutspec.utils import load_logger, basic_logger @@ -28,13 +28,13 @@ nucl_order5 = ['A', 'C', 'G', 'T', '-'] -def calc_phylocoefs(tree: PhyloTree): - tree_len = get_tree_len(tree, 'geom_mean') - logger.info(f'Tree len = {tree_len:.3f}') - phylocoefs = {tree.name: 1 - min(0.999, tree.get_closest_leaf()[1] / tree_len)} +def calc_phylocoefs(tree: Tree): + tree_height = get_tree_height(tree, 'geom_mean') + logger.info(f'Tree height = {tree_height:.3f}') + phylocoefs = {tree.name: 1 - min(0.999, tree.get_closest_leaf()[1] / tree_height)} for node in tree.iter_descendants(): _closest, d = node.get_closest_leaf() - phylocoefs[node.name] = 1 - min(0.99999, d / tree_len) + phylocoefs[node.name] = 1 - min(0.99999, d / tree_height) logger.info(f'Phylocoefs range: [{min(phylocoefs.values()):.3f}, {max(phylocoefs.values()):.3f}]') return phylocoefs @@ -342,7 +342,7 @@ def __init__( logger.info(f"Minimal probability for mutations to use: {proba_cutoff}") self.fp_format = np.float32 - self.tree = PhyloTree(path_to_tree, format=1) + self.tree = Tree(path_to_tree, format=1) logger.info( f"Tree loaded, number of leaf nodes: {len(self.tree)}, " f"total number of nodes: {len(self.tree.get_cached_content())}, " diff --git a/scripts/plot_spectra.py b/scripts/plot_spectra.py index 4ca8fe2..54c23f8 100644 --- a/scripts/plot_spectra.py +++ b/scripts/plot_spectra.py @@ -1,7 +1,5 @@ #!/usr/bin/env python3 -import os -import sys import re from functools import partial diff --git a/scripts/rename_internal_nodes.py b/scripts/rename_internal_nodes.py index d3be9ef..3557880 100644 --- a/scripts/rename_internal_nodes.py +++ b/scripts/rename_internal_nodes.py @@ -1,12 +1,12 @@ #!/usr/bin/env python3 import sys -from ete3 import PhyloTree, PhyloNode +from pymutspec.annotation.phylo_tree import Tree, TreeNode dist_formatter = "%0.8f" -def node_parent(node: PhyloNode): +def node_parent(node: TreeNode): try: return next(node.iter_ancestors()) except BaseException: @@ -20,8 +20,8 @@ def main(): except: print("ERROR\nUSAGE: script.py path_to_dist_tree path_to_named_tree path_to_out_tree", file=sys.stderr) - tree_dist = PhyloTree(path_to_dist_tree, format=0) - tree_named = PhyloTree(path_to_named_tree, format=8) + tree_dist = Tree(path_to_dist_tree, format=0) + tree_named = Tree(path_to_named_tree, format=8) nd = len(tree_dist.get_cached_content()) nn = len(tree_named.get_cached_content()) @@ -44,7 +44,7 @@ def main(): node_named = pa_named nwk = tree_dist.write(format=1, outfile=None, dist_formatter=dist_formatter) - nwk = nwk.replace(";", "ROOT;") # add ROOT label, that cannot be added with PhyloTree + nwk = nwk.replace(";", "ROOT;") # add ROOT label, that cannot be added with Tree with open(path_to_out, "w") as fout: fout.write(nwk) diff --git a/src/pymutspec/__init__.py b/src/pymutspec/__init__.py new file mode 100644 index 0000000..55581a9 --- /dev/null +++ b/src/pymutspec/__init__.py @@ -0,0 +1,8 @@ +from .draw import plot_mutspec, plot_mutspec12, plot_mutspec192 +from .annotation import ( + CodonAnnotation, mutations_summary, + calculate_mutspec, complete_sbs192_columns, + rev_comp, transcriptor +) + +__version__ = "0.0.15" \ No newline at end of file diff --git a/pymutspec/annotation/__init__.py b/src/pymutspec/annotation/__init__.py similarity index 86% rename from pymutspec/annotation/__init__.py rename to src/pymutspec/annotation/__init__.py index 3a9c792..be720b0 100644 --- a/pymutspec/annotation/__init__.py +++ b/src/pymutspec/annotation/__init__.py @@ -1,5 +1,6 @@ from .auxiliary import lbl2lbl_id, lbl_id2lbl, rev_comp, transcriptor -from .tree import iter_tree_edges, node_parent, calc_phylocoefs, get_ingroup_root, get_tree_len +from .phylo_tree import TreeNode, Tree +from .tree import iter_tree_edges, node_parent, calc_phylocoefs, get_ingroup_root, get_tree_height from .spectra import ( calculate_mutspec, collapse_mutspec, calc_edgewise_spectra, complete_sbs192_columns, jackknife_spectra_sampling, collapse_sbs192, diff --git a/src/pymutspec/annotation/auxiliary.py b/src/pymutspec/annotation/auxiliary.py new file mode 100644 index 0000000..4e641b7 --- /dev/null +++ b/src/pymutspec/annotation/auxiliary.py @@ -0,0 +1,86 @@ +transcriptor = str.maketrans("ACGT", "TGCA") + + +def rev_comp(mut: str): + """ + Return the reverse complement of a 192-component SBS mutation string. + + The input format is ``X[N>M]Y`` where ``X`` and ``Y`` are the flanking + nucleotides and ``N>M`` is the substitution. The function swaps the + flanking nucleotides and applies complement translation to all characters. + + Arguments + --------- + mut: str + SBS mutation string of the form ``X[N>M]Y``. + + Return + ------ + rev_comp_mut: str + Reverse-complemented mutation string. + """ + new_mut = mut[-1] + mut[1:-1] + mut[0] + new_mut = new_mut.translate(transcriptor) + return new_mut + + +def lbl_id2lbl(lbl_id: int) -> str: + """ + Convert a numeric label identifier to its string label. + + Arguments + --------- + lbl_id: int + Integer label code: 0 → ``'all'``, 1 → ``'syn'``, 2 → ``'ff'``. + + Return + ------ + lbl: str + Human-readable label string. + + Raises + ------ + NotImplementedError + If ``lbl_id`` is not 0, 1, or 2. + """ + if lbl_id == 0: + lbl = "all" + elif lbl_id == 1: + lbl = "syn" + elif lbl_id == 2: + lbl = "ff" + else: + raise NotImplementedError() + return lbl + + +def lbl2lbl_id(lbl: str) -> int: + """ + Convert a string label to its numeric identifier. + + Arguments + --------- + lbl: str + Label string; one of ``'all'``, ``'syn'``, ``'syn_c'``, + ``'ff'``, or ``'syn4f'``. + + Return + ------ + lbl_id: int + Integer label code: ``'all'`` → 0, ``'syn'``/``'syn_c'`` → 1, + ``'ff'``/``'syn4f'`` → 2. + + Raises + ------ + NotImplementedError + If ``lbl`` is not a recognised label string. + """ + if lbl == "all": + lbl_id = 0 + elif lbl == "syn" or lbl == "syn_c": + lbl_id = 1 + elif lbl == "ff" or lbl == "syn4f": + lbl_id = 2 + else: + raise NotImplementedError() + return lbl_id diff --git a/pymutspec/annotation/mut.py b/src/pymutspec/annotation/mut.py similarity index 98% rename from pymutspec/annotation/mut.py rename to src/pymutspec/annotation/mut.py index d6d7479..297ac96 100644 --- a/pymutspec/annotation/mut.py +++ b/src/pymutspec/annotation/mut.py @@ -1,5 +1,4 @@ import os -import sys from collections import defaultdict from typing import Set, Union, Dict, Iterable import multiprocessing as mp @@ -9,9 +8,9 @@ import pandas as pd from Bio.Data import CodonTable from Bio.Data.CodonTable import NCBICodonTableDNA -from ete3 import PhyloTree -from ..constants import * +from .phylo_tree import Tree +from ..constants import possible_nucls from ..utils import basic_logger from ..io import GenomeStatesTotal from .tree import iter_tree_edges, calc_phylocoefs @@ -439,9 +438,7 @@ def collect_exp_mut_freqs_proba( n = len(cds) if mask is not None and len(mask) != n: msg = f"Mask (len = {len(mask)}) must have same lenght as cds (len = {n})" - print(msg, file=sys.stderr) - # logger.error(msg) - # logger.info("Termination") + self.logger.error(msg) raise ValueError(msg) assert n % 3 == 0, "genomes length must be divisible by 3 (codon structure)" @@ -506,13 +503,17 @@ def collect_exp_muts_proba( n = len(cds) if mask is not None and len(mask) != n: msg = f"Mask (len = {len(mask)}) must have same lenght as cds (len = {n})" - print(msg, file=sys.stderr) - # logger.error(msg) - # logger.info("Termination") + self.logger.error(msg) raise ValueError(msg) - assert n % 3 == 0, "genomes length must be divisible by 3 (codon structure)" - assert 0 < phylocoef <= 1, "Evol coefficient must be between 0 and 1" + if n % 3 != 0: + self.logger.warning(f"genomes length ({n}) is not divisible by 3 (codon structure). Last codon will be skipped in syn, syn4f and pos3 modes") + n = n - (n % 3) + + if phylocoef <= 0 or phylocoef > 1: + msg = f"Evol coefficient must be between 0 and 1, but got {phylocoef}" + self.logger.error(msg) + raise ValueError(msg) labels = set(labels) data = [] @@ -975,7 +976,7 @@ def __init__( self.logger.info(f"Minimal probability for mutations to use: {proba_cutoff}") self.fp_format = np.float32 - self.tree = PhyloTree(path_to_tree, format=1) + self.tree = Tree(path_to_tree, format=1) self.logger.info( f"Tree loaded, number of leaf nodes: {len(self.tree)}, " f"total number of nodes: {len(self.tree.get_cached_content())}, " @@ -1097,7 +1098,7 @@ def mutations_summary(mutations: pd.DataFrame, gene_col=None, proba_col=None, ge --------- mutations: pd.DataFrame table must contain at least 2 columns: - - Mut: str; Pattern: '[ACGT]\[[ACGT]>[ACGT]\][ACGT]' + - Mut: str; Pattern: ``[ACGT]\\[[ACGT]>[ACGT]\\][ACGT]`` - Label: int; [-3, 2]. See CodonAnnotation.get_mut_type - $gene_col, optional. If gene_col=None annotation will be formed on full mutations without genes splitting - $proba_col, optional. If proba_col=None each row of table assumed to be mutation, else probabilities will be used @@ -1115,7 +1116,7 @@ def mutations_summary(mutations: pd.DataFrame, gene_col=None, proba_col=None, ge table with mutations annotation """ mutations = mutations.copy() - mut_pattern = "[ACGT]\[[ACGT]>[ACGT]\][ACGT]" + mut_pattern = r"[ACGT]\[[ACGT]>[ACGT]\][ACGT]" label_mapper = { -3: "6Stop to stop", -2: "4Stop loss", diff --git a/src/pymutspec/annotation/phylo_tree.py b/src/pymutspec/annotation/phylo_tree.py new file mode 100644 index 0000000..8cfbae9 --- /dev/null +++ b/src/pymutspec/annotation/phylo_tree.py @@ -0,0 +1,220 @@ +""" +Custom phylogenetic tree classes replacing ete3 dependency. +Uses BioPython for newick format parsing. +""" +from io import StringIO + +from Bio import Phylo as _BioPhylo + + +class TreeNode: + """ + A phylogenetic tree node with an ete3-compatible interface. + Instances form a tree through parent/child references. + """ + + def __init__(self, name="", dist=0.0): + self.name = name if name is not None else "" + self.dist = dist if dist is not None else 0.0 + self.children = [] + self._parent = None + + # ------------------------------------------------------------------ + # Tree navigation + # ------------------------------------------------------------------ + + def is_leaf(self): + return len(self.children) == 0 + + def traverse(self): + """Yield all nodes (self first, then descendants).""" + yield self + for child in self.children: + yield from child.traverse() + + def iter_descendants(self): + """Yield all descendant nodes (not self).""" + for child in self.children: + yield child + yield from child.iter_descendants() + + def iter_ancestors(self): + """Yield ancestor nodes from parent to root.""" + node = self._parent + while node is not None: + yield node + node = node._parent + + def iter_leaves(self): + """Yield all leaf nodes reachable from self.""" + if self.is_leaf(): + yield self + else: + for child in self.children: + yield from child.iter_leaves() + + # ------------------------------------------------------------------ + # Distance computation + # ------------------------------------------------------------------ + + def get_distance(self, target): + """Return cumulative branch length from self to *target* (descendant).""" + def _find(node, target, acc): + if node is target: + return acc + for child in node.children: + result = _find(child, target, acc + child.dist) + if result is not None: + return result + return None + + d = _find(self, target, 0.0) + if d is None: + raise ValueError(f"Node '{target.name}' is not a descendant of '{self.name}'") + return d + + def get_farthest_leaf(self): + """Return *(leaf, distance)* for the leaf farthest from self.""" + best_leaf, best_dist = None, -1.0 + for leaf in self.iter_leaves(): + d = self.get_distance(leaf) + if d > best_dist: + best_dist = d + best_leaf = leaf + return best_leaf, best_dist + + def get_closest_leaf(self): + """Return *(leaf, distance)* for the leaf closest to self.""" + best_leaf, best_dist = None, float("inf") + for leaf in self.iter_leaves(): + d = self.get_distance(leaf) + if d < best_dist: + best_dist = d + best_leaf = leaf + return best_leaf, best_dist + + # ------------------------------------------------------------------ + # Node lookup + # ------------------------------------------------------------------ + + def get_cached_content(self): + """Return a dict keyed by every node reachable from self (ete3 compat).""" + return {node: None for node in self.traverse()} + + def search_nodes(self, name=None, **kwargs): + """Return list of all nodes whose attributes match the given criteria.""" + results = [] + for node in self.traverse(): + if name is not None and node.name != name: + continue + if all(getattr(node, k, None) == v for k, v in kwargs.items()): + results.append(node) + return results + + def iter_search_nodes(self, name=None, **kwargs): + """Yield all nodes whose attributes match the given criteria.""" + for node in self.traverse(): + if name is not None and node.name != name: + continue + if all(getattr(node, k, None) == v for k, v in kwargs.items()): + yield node + + # ------------------------------------------------------------------ + # Newick serialisation + # ------------------------------------------------------------------ + + def _to_newick(self, dist_formatter=None): + fmt = dist_formatter if dist_formatter else "%g" + if self.is_leaf(): + return f"{self.name}:{fmt % self.dist}" + children_str = ",".join(c._to_newick(dist_formatter) for c in self.children) + return f"({children_str}){self.name}:{fmt % self.dist}" + + def write(self, format=1, outfile=None, dist_formatter=None): # noqa: A002 + """ + Serialise tree to newick. + + Parameters + ---------- + format : int + Newick format hint (kept for API compatibility with ete3; currently + ignored – the tree is always written in a standard newick format + that includes all node names and branch lengths). + outfile : str or None + If given, also write the string to this file path. + dist_formatter : str or None + printf-style format string for branch lengths (e.g. ``"%0.8f"``). + + Returns + ------- + str + Newick string. + """ + children_str = ",".join(c._to_newick(dist_formatter) for c in self.children) + fmt = dist_formatter if dist_formatter else "%g" + nwk = f"({children_str}){self.name}:{fmt % self.dist};" + if outfile: + with open(outfile, "w") as fh: + fh.write(nwk) + return nwk + + def __len__(self): + """Return the number of leaf nodes (mirrors ete3 behaviour).""" + return sum(1 for _ in self.iter_leaves()) + + def __repr__(self): + return f"TreeNode(name={self.name!r}, dist={self.dist}, children={len(self.children)})" + + +# --------------------------------------------------------------------------- +# Internal helpers +# --------------------------------------------------------------------------- + +def _bio_clade_to_node(clade): + """Recursively convert a BioPython Clade into a TreeNode tree.""" + node = TreeNode( + name=clade.name if clade.name is not None else "", + dist=clade.branch_length if clade.branch_length is not None else 0.0, + ) + for child_clade in clade.clades: + child_node = _bio_clade_to_node(child_clade) + child_node._parent = node + node.children.append(child_node) + return node + + +# --------------------------------------------------------------------------- +# Public constructor – mirrors ete3's PhyloTree(path, format=N) +# --------------------------------------------------------------------------- + +class Tree(TreeNode): + """ + Load a phylogenetic tree from a newick file. + + Parameters + ---------- + newick_path : str + Path to a newick-format tree file. + format : int + Newick format hint (kept for API compatibility with ete3; this + parameter is currently ignored – BioPython's newick parser handles + all common newick variants automatically). + """ + + def __init__(self, newick_path, format=1): # noqa: A002 + with open(newick_path) as fh: + tree_str = fh.read().strip() + + bio_tree = _BioPhylo.read(StringIO(tree_str), "newick") + root = _bio_clade_to_node(bio_tree.root) + + # Initialise self as the root node (TreeNode.__init__ not called via + # super because we copy attributes from the parsed root directly). + self.name = root.name + self.dist = root.dist + self.children = root.children + self._parent = None + # Re-point children's _parent to self (they currently point to root). + for child in self.children: + child._parent = self + diff --git a/pymutspec/annotation/spectra.py b/src/pymutspec/annotation/spectra.py similarity index 62% rename from pymutspec/annotation/spectra.py rename to src/pymutspec/annotation/spectra.py index 5d2e55c..853dcc0 100644 --- a/pymutspec/annotation/spectra.py +++ b/src/pymutspec/annotation/spectra.py @@ -1,5 +1,5 @@ from sys import stderr -from typing import Set, Union, Dict, Iterable +from typing import Dict import numpy as np import pandas as pd @@ -30,7 +30,7 @@ def calculate_mutspec( --------- obs_muts: pd.DataFrame table containing mutations with annotation; table must contain 2 columns: - - Mut: str; Pattern: '[ACGT]\[[ACGT]>[ACGT]\][ACGT]' + - Mut: str; Pattern: ``[ACGT]\\[[ACGT]>[ACGT]\\][ACGT]`` - ProbaFull (optional, only for use_proba=True) - probability of mutation exp_muts: dict[str, float] @@ -145,6 +145,25 @@ def get_iqr_bounds(series: pd.Series): def filter_outlier_branches(obs_df: pd.DataFrame, use_proba=True): + """ + Remove branches with an outlier-high number of observed mutations. + + Outliers are identified using the IQR method: branches whose mutation + count exceeds ``Q3 + 1.5 * IQR`` are dropped. + + Arguments + --------- + obs_df: pd.DataFrame + Observed-mutations table containing at least the columns + ``'AltNode'``, ``'Mut'``, and optionally ``'ProbaMut'``. + use_proba: bool + If ``True`` sum ``'ProbaMut'`` per branch; otherwise count rows. + + Return + ------ + obs_df_flt: pd.DataFrame + Filtered mutations table with outlier branches removed. + """ if use_proba: edge_nobs = obs_df.groupby('AltNode')['ProbaMut'].sum() else: @@ -158,9 +177,35 @@ def filter_outlier_branches(obs_df: pd.DataFrame, use_proba=True): def collapse_mutspec(ms192: pd.DataFrame): - assert ms192.shape[0] == 192 + """ + Collapse a 192-component spectrum to 96 components using reverse complement. + + Mutations on the ``A``/``G`` strand are reverse-complemented so that all + substitutions are expressed relative to the pyrimidine base (``C`` or + ``T``), then the ``ObsFr`` and ``ExpFr`` columns are summed for matching + contexts, yielding a 96-component spectrum. + + Arguments + --------- + ms192: pd.DataFrame + 192-component spectrum table. Must contain columns ``'Mut'``, + ``'ObsFr'``, and ``'ExpFr'``, and must have exactly 192 rows. + + Return + ------ + ms96: pd.DataFrame + 96-component spectrum with columns ``'ObsFr'``, ``'ExpFr'``, + ``'RawMutSpec'``, and ``'MutSpec'`` (normalised to sum to 1). + + Raises + ------ + AssertionError + If ``ms192`` does not have exactly 192 rows or is missing required + columns. + """ + assert ms192.shape[0] == 192, f"Expected 192 rows, got {ms192.shape[0]}" for c in ["Mut", "ObsFr", "ExpFr"]: - assert c in ms192.columns + assert c in ms192.columns, f"Required column '{c}' not found in ms192" ms1 = ms192[ms192["Mut"].str.get(2).isin(list("CT"))] ms2 = ms192[ms192["Mut"].str.get(2).isin(list("AG"))] @@ -174,6 +219,23 @@ def collapse_mutspec(ms192: pd.DataFrame): def complete_sbs192_columns(df: pd.DataFrame): + """ + Ensure a DataFrame has all 192 SBS columns, filling missing ones with 0. + + The resulting DataFrame is reordered so its columns follow the canonical + ``possible_sbs192`` order. + + Arguments + --------- + df: pd.DataFrame + DataFrame whose columns are a (possibly incomplete) subset of the 192 + SBS mutation types. + + Return + ------ + df: pd.DataFrame + DataFrame with exactly 192 columns in canonical order. + """ df = df.copy() if len(df.columns) != 192: for sbs192 in possible_sbs192_set.difference(df.columns.values): @@ -183,7 +245,36 @@ def complete_sbs192_columns(df: pd.DataFrame): def collapse_sbs192(df: pd.DataFrame, to=12): - assert (df.columns == possible_sbs192).all() + """ + Sum a 192-component SBS DataFrame into a 12-component representation. + + Each 192-component mutation type is mapped to its 12-component base + substitution (the middle three characters, e.g. ``'C>A'``), and the + values are accumulated. + + Arguments + --------- + df: pd.DataFrame + DataFrame with columns equal to ``possible_sbs192`` in canonical order. + Each row typically represents one sample or branch. + to: int + Target number of components. Currently only ``12`` is supported. + + Return + ------ + df12: pd.DataFrame + DataFrame with 12 columns corresponding to the 12 base substitution + types in ``possible_sbs12`` order. + + Raises + ------ + AssertionError + If ``df.columns`` does not match ``possible_sbs192``. + NotImplementedError + If ``to`` is not ``12``. + """ + assert (df.columns == possible_sbs192).all(), \ + "DataFrame columns must match possible_sbs192 in canonical order" df = df.copy() if to == 12: for sbs192 in possible_sbs192: @@ -199,6 +290,36 @@ def collapse_sbs192(df: pd.DataFrame, to=12): def jackknife_spectra_sampling(obs: pd.DataFrame, exp: pd.DataFrame, frac=0.5, n=1000): + """ + Estimate spectrum variability via jackknife resampling of tree branches. + + On each iteration a random subset of branches (edges) is drawn without + replacement and a per-branch spectrum ratio ``obs / exp`` is computed. + The resulting collection of spectra can be used to derive confidence + intervals. + + Arguments + --------- + obs: pd.DataFrame + Observed mutations. Either a pre-pivoted wide DataFrame with 192 + SBS columns and a ``(RefNode, AltNode)`` MultiIndex, or a long-format + DataFrame with columns ``'AltNode'``, ``'RefNode'``, ``'Mut'``, + and ``'ProbaFull'``. + exp: pd.DataFrame + Expected mutation frequencies. Either a pre-pivoted wide DataFrame + with 192 SBS columns and a ``Node`` index, or a long-format DataFrame + with columns ``'Node'``, ``'Mut'``, and ``'Proba'``. + frac: float + Fraction of branches to sample on each iteration. + n: int + Number of jackknife iterations. + + Return + ------ + spectra: pd.DataFrame + DataFrame of shape ``(n, 192)`` where each row is the spectrum + computed from one jackknife sample. + """ if len(obs.columns) == 192 and \ (obs.columns == possible_sbs192).all() and \ (exp.columns == possible_sbs192).all(): @@ -235,28 +356,51 @@ def jackknife_spectra_sampling(obs: pd.DataFrame, exp: pd.DataFrame, frac=0.5, n return pd.DataFrame(spectra).fillna(0.) -def collapse_sbs192(df: pd.DataFrame, to=12): - assert (df.columns == possible_sbs192).all() - df = df.copy() - if to == 12: - for sbs192 in possible_sbs192: - sbs12 = sbs192[2:5] - if sbs12 in df.columns.values: - df[sbs12] += df[sbs192] - else: - df[sbs12] = df[sbs192] - - return df[possible_sbs12] - else: - raise NotImplementedError() - - def calc_edgewise_spectra( obs: pd.DataFrame, exp: pd.DataFrame, nmtypes_cutoff=10, nobs_cuttof=10, collapse_to_12=False, scale=True, both_12_and_192=False ): + """ + Calculate per-branch (edge-wise) mutational spectra. + + For each tree branch the observed mutation counts are divided by the + expected frequencies of the reference (parent) node, yielding a + branch-specific spectrum. + + Arguments + --------- + obs: pd.DataFrame + Observed mutations. Either a pre-pivoted wide DataFrame with 192 + SBS columns and a ``(RefNode, AltNode)`` MultiIndex, or a long-format + DataFrame with columns ``'RefNode'``, ``'AltNode'``, ``'Mut'``, and + ``'ProbaFull'``. + exp: pd.DataFrame + Expected mutation frequencies. Either a pre-pivoted wide DataFrame + with 192 SBS columns and a ``Node`` index, or a long-format DataFrame + with columns ``'Node'``, ``'Mut'``, and ``'Proba'``. + nmtypes_cutoff: int + Minimum number of distinct mutation types a branch must have to be + retained (only applied when ``collapse_to_12=False``). + nobs_cuttof: int + Minimum total observed mutations a branch must have to be retained + (only applied when ``collapse_to_12=False``). + collapse_to_12: bool + If ``True`` collapse the 192-component spectra to 12 components before + returning. + scale: bool + If ``True`` normalise each branch spectrum to sum to 1. + both_12_and_192: bool + If ``True`` return both 12- and 192-component spectra as a tuple + ``(spectra12, spectra192)``. + + Return + ------ + spectra: pd.DataFrame or tuple[pd.DataFrame, pd.DataFrame] + Branch-wise spectrum DataFrame (or tuple of two DataFrames when + ``both_12_and_192=True``). + """ if len(obs.columns) == 192 and \ (obs.columns == possible_sbs192).all() and \ (exp.columns == possible_sbs192).all(): @@ -321,7 +465,24 @@ def calc_edgewise_spectra( def get_cossim(a: pd.DataFrame, b: pd.DataFrame): - assert (a.columns == b.columns).all() + """ + Compute row-wise cosine similarity between two aligned DataFrames. + + Only rows present in both DataFrames (intersection of indices) are used. + + Arguments + --------- + a: pd.DataFrame + First DataFrame; columns must match those of *b*. + b: pd.DataFrame + Second DataFrame; columns must match those of *a*. + + Return + ------ + cossim: pd.Series + Cosine similarity for each shared index, ranging from -1 to 1. + Returns an empty Series if the indices do not overlap. + """ common_index = a.index.intersection(b.index) if len(common_index) == 0: @@ -338,7 +499,24 @@ def get_cossim(a: pd.DataFrame, b: pd.DataFrame): def get_eucdist(a: pd.DataFrame, b: pd.DataFrame): - assert (a.columns == b.columns).all() + """ + Compute row-wise Euclidean distance between two aligned DataFrames. + + Only rows present in both DataFrames (intersection of indices) are used. + + Arguments + --------- + a: pd.DataFrame + First DataFrame; columns must match those of *b*. + b: pd.DataFrame + Second DataFrame; columns must match those of *a*. + + Return + ------ + d: pd.Series + Euclidean distance for each shared index. + Returns an empty Series if the indices do not overlap. + """ common_index = a.index.intersection(b.index) if len(common_index) == 0: diff --git a/src/pymutspec/annotation/tree.py b/src/pymutspec/annotation/tree.py new file mode 100644 index 0000000..382ce8a --- /dev/null +++ b/src/pymutspec/annotation/tree.py @@ -0,0 +1,179 @@ +from queue import Queue +from statistics import geometric_mean + +import numpy as np + +# TODO merge with phylo_tree.py and remove this file; fix tests after these changes + +def node_parent(node): + """ + TODO remove this function and use node._parent directly AND rename node._parent to node.parent + + Return the parent node of *node*, or ``None`` if *node* is the root. + + Arguments + --------- + node + A node in a phylogenetic tree. + + Return + ------ + parent or None + The immediate ancestor of *node*, or ``None`` when *node* has no + ancestors (i.e. it is the root). + """ + try: + return next(node.iter_ancestors()) + except BaseException: + return None + + +def iter_tree_edges(tree): + """ + TODO integrate to Tree class + + Iterate over all directed edges (parent → child) in the tree via BFS. + + The root node itself is skipped; every other node produces exactly one + ``(ref_node, alt_node)`` pair where *ref_node* is the parent and + *alt_node* is the child. + + Arguments + --------- + tree + Rooted phylogenetic tree. + + Yields + ------ + ref_node + Parent (reference) node of the edge. + alt_node + Child (alternative) node of the edge. + """ + discovered_nodes = set() + discovered_nodes.add(tree.name) + Q = Queue() + Q.put(tree) + + while not Q.empty(): + cur_node = Q.get() + for child in cur_node.children: + Q.put(child) + + if cur_node.name not in discovered_nodes: + discovered_nodes.add(cur_node.name) + alt_node = cur_node + ref_node = node_parent(alt_node) + yield ref_node, alt_node + + +def get_tree_height(tree, mode='geom_mean'): + """ + Return the characteristic length of a (sub)tree as the distance from + the root to its leaves. + + Arguments + --------- + tree + Rooted phylogenetic tree or subtree. Must not be named ``'ROOT'``. + mode: str + Aggregation method over leaf distances. One of: + + - ``'mean'`` – arithmetic mean of leaf distances + - ``'geom_mean'`` – geometric mean of leaf distances (default) + - ``'max'`` – distance to the farthest leaf + + Return + ------ + tree_len: float + Characteristic length of the tree. + + Raises + ------ + TypeError + If *mode* is not one of the accepted values. + """ + if mode == 'max': + _, md = tree.get_farthest_leaf() + elif mode in ['mean', 'geom_mean']: + distances_to_leaves = [] + for leaf in tree.iter_leaves(): + d = tree.get_distance(leaf) + distances_to_leaves.append(d) + + if mode == 'mean': + md = np.mean(distances_to_leaves) + elif mode == 'geom_mean': + md = geometric_mean(distances_to_leaves) + + else: + raise TypeError("mode must be 'mean', 'geom_mean' or 'max'") + + return md + + +def get_ingroup_root(tree): + """ + Return the ingroup root of a binary rooted tree that contains an outgroup. + + The function assumes the tree root has exactly two children, one of which + is a leaf (the outgroup). If no leaf child is found the tree root itself + is returned. + + Arguments + --------- + tree + Rooted binary tree with an outgroup leaf attached to the root. + + Return + ------ + ingrp + Root of the ingroup clade. + + Raises + ------ + AssertionError + If the tree root does not have exactly two children. + """ + assert len(tree.children) == 2, 'Tree must be binary' + found_outgroup = False + for node in tree.children: + if node.is_leaf(): + found_outgroup = True + else: + ingrp = node + + if found_outgroup: + return ingrp + else: + return tree + + +def calc_phylocoefs(tree): + """ + Calculate a phylogenetic coefficient for every node in the tree. + + The coefficient for a node is ``1 - d / tree_len``, where *d* is the + distance from the node to its closest leaf and *tree_len* is the + geometric-mean leaf distance of the ingroup root. Values are capped so + that the minimum coefficient is > 0 (i.e. ``d / tree_len`` is capped + at 0.99999). + + Arguments + --------- + tree + Rooted binary phylogenetic tree (with an outgroup leaf, optional). + + Return + ------ + phylocoefs: dict[str, float] + Mapping from node name to its phylogenetic coefficient. + """ + ingroup = get_ingroup_root(tree) + tree_height = get_tree_height(ingroup, 'geom_mean') + root_phylocoef = 1 - min(0.999, ingroup.get_closest_leaf()[1] / tree_height) + phylocoefs = {ingroup.name: root_phylocoef} + for node in ingroup.iter_descendants(): + _closest, d = node.get_closest_leaf() + phylocoefs[node.name] = 1 - min(0.99999, d / tree_height) + return phylocoefs diff --git a/pymutspec/constants/__init__.py b/src/pymutspec/constants/__init__.py similarity index 100% rename from pymutspec/constants/__init__.py rename to src/pymutspec/constants/__init__.py diff --git a/pymutspec/constants/sbs.py b/src/pymutspec/constants/sbs.py similarity index 100% rename from pymutspec/constants/sbs.py rename to src/pymutspec/constants/sbs.py diff --git a/src/pymutspec/draw/__init__.py b/src/pymutspec/draw/__init__.py new file mode 100644 index 0000000..d6b0689 --- /dev/null +++ b/src/pymutspec/draw/__init__.py @@ -0,0 +1,4 @@ +from .spectra import plot_mutspec, plot_mutspec12, plot_mutspec192 + +def plot_mutspec192box(*args, **kwargs): + print("WARNING: the function is removed! Use plot_mutspec192(style='box') instead") diff --git a/src/pymutspec/draw/sbs_orders.py b/src/pymutspec/draw/sbs_orders.py new file mode 100644 index 0000000..9aad453 --- /dev/null +++ b/src/pymutspec/draw/sbs_orders.py @@ -0,0 +1,22 @@ +# import pandas as pd + +# from ..constants import possible_sbs192 +# from ..annotation import rev_comp, transcriptor + +# kk_lbls = "A>C A>G A>T C>T G>C G>T".split() +# cosmic_lbls = "C>A C>G C>T T>A T>C T>G".split() + +# df = pd.DataFrame({"sbs": possible_sbs192}) +# df["sbs_base"] = df["sbs"].str.slice(2, 5) +# df["sbs_base_revcomp"] = df["sbs_base"].str.translate(transcriptor) +# df["sbs_revcomp"] = df["sbs"].apply(rev_comp) +# df["is_cosmic"] = df["sbs_base"].isin(cosmic_lbls) +# df["is_kk"] = df["sbs_base"].isin(kk_lbls) +# df["sbs_base_for_sorting_kp"] = df.apply( +# lambda x: x.sbs_base + "1" if x.is_cosmic else x.sbs_base_revcomp + "2", axis=1) +# df["sbs_for_ordering_kk"] = df.apply(lambda x: x.sbs if x.is_kk else x.sbs_revcomp, axis=1) +# df["sbs_for_ordering_kp"] = df.apply(lambda x: x.sbs if x.is_cosmic else x.sbs_revcomp, axis=1) + +# ordered_sbs192_kp = list(df.sort_values(["sbs_base_for_sorting_kp", "sbs_for_ordering_kp"]).sbs.values) +# ordered_sbs192_kk = list(df.sort_values(["sbs_base", "sbs_for_ordering_kk"]).sbs.values) +# del df diff --git a/src/pymutspec/draw/spectra.py b/src/pymutspec/draw/spectra.py new file mode 100644 index 0000000..9758b4d --- /dev/null +++ b/src/pymutspec/draw/spectra.py @@ -0,0 +1,515 @@ +""" +Functionality to plot mutational spectrums +""" + +from typing import Iterable + +import pandas as pd +import matplotlib.pyplot as plt +import seaborn as sns + +from ..constants import possible_sbs192 + +ordered_sbs12 = ["C>A", "G>T", "C>G", "G>C", "C>T", "G>A", + "T>A", "A>T", "T>C", "A>G", "T>G", "A>C"] +ordered_sbs192 = [ + 'A[C>A]A', 'A[C>A]C', 'A[C>A]G', 'A[C>A]T', 'C[C>A]A', 'C[C>A]C', + 'C[C>A]G', 'C[C>A]T', 'G[C>A]A', 'G[C>A]C', 'G[C>A]G', 'G[C>A]T', + 'T[C>A]A', 'T[C>A]C', 'T[C>A]G', 'T[C>A]T', 'T[G>T]T', 'G[G>T]T', + 'C[G>T]T', 'A[G>T]T', 'T[G>T]G', 'G[G>T]G', 'C[G>T]G', 'A[G>T]G', + 'T[G>T]C', 'G[G>T]C', 'C[G>T]C', 'A[G>T]C', 'T[G>T]A', 'G[G>T]A', + 'C[G>T]A', 'A[G>T]A', 'A[C>G]A', 'A[C>G]C', 'A[C>G]G', 'A[C>G]T', + 'C[C>G]A', 'C[C>G]C', 'C[C>G]G', 'C[C>G]T', 'G[C>G]A', 'G[C>G]C', + 'G[C>G]G', 'G[C>G]T', 'T[C>G]A', 'T[C>G]C', 'T[C>G]G', 'T[C>G]T', + 'T[G>C]T', 'G[G>C]T', 'C[G>C]T', 'A[G>C]T', 'T[G>C]G', 'G[G>C]G', + 'C[G>C]G', 'A[G>C]G', 'T[G>C]C', 'G[G>C]C', 'C[G>C]C', 'A[G>C]C', + 'T[G>C]A', 'G[G>C]A', 'C[G>C]A', 'A[G>C]A', 'A[C>T]A', 'A[C>T]C', + 'A[C>T]G', 'A[C>T]T', 'C[C>T]A', 'C[C>T]C', 'C[C>T]G', 'C[C>T]T', + 'G[C>T]A', 'G[C>T]C', 'G[C>T]G', 'G[C>T]T', 'T[C>T]A', 'T[C>T]C', + 'T[C>T]G', 'T[C>T]T', 'T[G>A]T', 'G[G>A]T', 'C[G>A]T', 'A[G>A]T', + 'T[G>A]G', 'G[G>A]G', 'C[G>A]G', 'A[G>A]G', 'T[G>A]C', 'G[G>A]C', + 'C[G>A]C', 'A[G>A]C', 'T[G>A]A', 'G[G>A]A', 'C[G>A]A', 'A[G>A]A', + 'A[T>A]A', 'A[T>A]C', 'A[T>A]G', 'A[T>A]T', 'C[T>A]A', 'C[T>A]C', + 'C[T>A]G', 'C[T>A]T', 'G[T>A]A', 'G[T>A]C', 'G[T>A]G', 'G[T>A]T', + 'T[T>A]A', 'T[T>A]C', 'T[T>A]G', 'T[T>A]T', 'T[A>T]T', 'G[A>T]T', + 'C[A>T]T', 'A[A>T]T', 'T[A>T]G', 'G[A>T]G', 'C[A>T]G', 'A[A>T]G', + 'T[A>T]C', 'G[A>T]C', 'C[A>T]C', 'A[A>T]C', 'T[A>T]A', 'G[A>T]A', + 'C[A>T]A', 'A[A>T]A', 'A[T>C]A', 'A[T>C]C', 'A[T>C]G', 'A[T>C]T', + 'C[T>C]A', 'C[T>C]C', 'C[T>C]G', 'C[T>C]T', 'G[T>C]A', 'G[T>C]C', + 'G[T>C]G', 'G[T>C]T', 'T[T>C]A', 'T[T>C]C', 'T[T>C]G', 'T[T>C]T', + 'T[A>G]T', 'G[A>G]T', 'C[A>G]T', 'A[A>G]T', 'T[A>G]G', 'G[A>G]G', + 'C[A>G]G', 'A[A>G]G', 'T[A>G]C', 'G[A>G]C', 'C[A>G]C', 'A[A>G]C', + 'T[A>G]A', 'G[A>G]A', 'C[A>G]A', 'A[A>G]A', 'A[T>G]A', 'A[T>G]C', + 'A[T>G]G', 'A[T>G]T', 'C[T>G]A', 'C[T>G]C', 'C[T>G]G', 'C[T>G]T', + 'G[T>G]A', 'G[T>G]C', 'G[T>G]G', 'G[T>G]T', 'T[T>G]A', 'T[T>G]C', + 'T[T>G]G', 'T[T>G]T', 'T[A>C]T', 'G[A>C]T', 'C[A>C]T', 'A[A>C]T', + 'T[A>C]G', 'G[A>C]G', 'C[A>C]G', 'A[A>C]G', 'T[A>C]C', 'G[A>C]C', + 'C[A>C]C', 'A[A>C]C', 'T[A>C]A', 'G[A>C]A', 'C[A>C]A', 'A[A>C]A' +] +ordered_sbs192_kp = ordered_sbs192 + +# KK-style ordering: substitutions are grouped by base type and sorted using +# the SBS itself for kk_lbls, or its reverse-complement otherwise. +_kk_lbl_set = set("A>C A>G A>T C>T G>C G>T".split()) +_transcriptor = str.maketrans("ACGT", "TGCA") + + +def _sbs192_rev_comp(sbs: str) -> str: + """Return the reverse complement of a 192-component SBS string. + + The input must be a 7-character string of the form ``X[N>M]Y`` where + ``X`` and ``Y`` are single flanking nucleotides and ``N>M`` is the + substitution (e.g. ``'A[C>A]T'``). + """ + return (sbs[-1] + sbs[1:-1] + sbs[0]).translate(_transcriptor) + + +ordered_sbs192_kk = sorted( + possible_sbs192, + key=lambda sbs: (sbs[2:5], sbs if sbs[2:5] in _kk_lbl_set else _sbs192_rev_comp(sbs)), +) + +color_mapping12 = { + "C>A": "deepskyblue", + "G>T": "deepskyblue", + "C>G": "black", + "G>C": "black", + "C>T": "red", + "G>A": "red", + "T>A": "silver", + "A>T": "silver", + "T>C": "yellowgreen", + "A>G": "yellowgreen", + "T>G": "pink", + "A>C": "pink", +} +color_mapping192 = {} +for _sbs192 in ordered_sbs192: + _sbs12 = _sbs192[2:5] + color_mapping192[_sbs192] = color_mapping12[_sbs12] + + +def _prepare_nice_labels(sbs192: Iterable[str], kk=False): + """ + Build a display-friendly label list for 192-component SBS axes. + + A short separator string (underscores) is inserted between groups of + substitutions that share the same base substitution type. + + Arguments + --------- + sbs192: iterable of str + Ordered list of 192-component SBS strings (e.g. ``'A[C>A]C'``). + kk: bool + If ``True``, use the compact KK-style label format + ``'CA: ACA'`` instead of the full COSMIC string. + + Return + ------ + labels: list of str + Label strings suitable for use as tick labels, with separator + entries inserted between substitution groups. + """ + _nice_sbs = [] + prev = None + for i, sbs in enumerate(sbs192, 1): + if prev is not None and sbs[2:5] != prev[2:5]: + _nice_sbs.append("_" * (i // 10)) + sbs_nice = sbs[2] + sbs[4] + ": " + sbs[0] + sbs[2] + sbs[-1] if kk else sbs + _nice_sbs.append(sbs_nice) + prev = sbs + return _nice_sbs + + +def plot_mutspec( + mutspec: pd.DataFrame, + spectra_col="MutSpec", + title="Spectrum", + ylabel=None, + figsize=None, + style="bar", + sbs_kind=12, + sbs_order=None, + labels_style="cosmic", + savepath=None, + fontname=None, + ticksize=8, + titlesize=14, + ylabelsize=12, + bar_width=0.8, + show=True, + dpi=300, + ax=None, + **kwargs, + ): + """ + General plotting function for mutational spectra supporting + 12- and 192-component spectra. + + Arguments + --------- + mutspec: pd.DataFrame + Table containing at least a ``'Mut'`` column and a column named + *spectra_col* with spectrum values. + spectra_col: str + Column name used for the y-axis values. + title: str + Plot title. + ylabel: str or None + Y-axis label. + figsize: tuple or None + Figure size ``(width, height)`` in inches. Defaults to ``(24, 8)`` + for 192-component and ``(6, 4)`` for 12-component spectra. + style: str + ``'bar'`` for a bar plot or ``'box'`` for a box plot. + sbs_kind: int + Number of SBS components: ``12`` or ``192``. + sbs_order: list or None + Custom order for SBS labels on the x-axis. + labels_style: str + Style for 192-component tick labels. One of: + + - ``'cosmic'`` – full COSMIC-format strings, e.g. ``'A[C>A]T'`` + (default) + - ``'long'`` – same as ``'cosmic'`` + - ``'kk'`` – compact KK-style labels, e.g. ``'CA: ACT'`` + savepath: str or None + File path to save the figure. If ``None`` the figure is not saved. + fontname: str or None + Font family for all text elements. + ticksize: int + Font size for tick labels. + titlesize: int + Font size for the title. + ylabelsize: int + Font size for the y-axis label. + bar_width: float + Width of bars in bar plots. + show: bool + If ``True`` call ``plt.show()``; otherwise close the figure. + dpi: int + Resolution for saved figures. + ax: matplotlib.axes.Axes or None + Axes to draw on. A new figure is created when ``None``. + **kwargs + Additional keyword arguments forwarded to the underlying seaborn + plot function. + + Return + ------ + ax: matplotlib.axes.Axes + The axes containing the plot. + """ + if "filepath" in kwargs: + savepath = kwargs.pop("filepath") + + is_192 = (sbs_kind == 192) + kk_labels = (labels_style == "kk") + + # Defaults + if figsize is None: + figsize = (24, 8) if is_192 else (6, 4) + + ms = mutspec.copy() + + if is_192: + sbs_order = sbs_order or ordered_sbs192 + order = _prepare_nice_labels(sbs_order, kk=kk_labels) + if kk_labels: + # Map COSMIC SBS strings to KK-style display labels and rebuild palette + sbs_to_kk = { + sbs: sbs[2] + sbs[4] + ": " + sbs[0] + sbs[2] + sbs[-1] + for sbs in possible_sbs192 + } + ms['Mut'] = ms['Mut'].map(sbs_to_kk) + palette = {sbs_to_kk[sbs]: color_mapping192[sbs] for sbs in possible_sbs192} + else: + palette = color_mapping192 + tick_rotation = 90 + else: + order = sbs_order or ordered_sbs12 + palette = color_mapping12 + tick_rotation = 0 + + if ax is None: + fig = plt.figure(figsize=figsize) + ax = fig.gca() + + if style == "bar": + _cols = set(ms.columns) + if 'MutSpec_median' in _cols and 'MutSpec_q05' in _cols and 'MutSpec_q95' in _cols: + sns.barplot( + ms, x='Mut', y='MutSpec', hue='Mut', legend=False, width=bar_width, + order=order, palette=palette, err_kws={'linewidth': 1}, ax=ax, **kwargs) + + mutspec_index = ms.set_index('Mut') + + # build y and yerr aligned with order; support '_' separators in 192 order + ymed, yerr_min, yerr_max = [], [], [] + for mt in order: + if isinstance(mt, str) and mt.startswith('_'): + ymed.append(0.) + yerr_min.append(0.) + yerr_max.append(0.) + else: + row = mutspec_index.loc[mt] + ymed.append(row['MutSpec_median']) + # convert to relative error for plotting + yerr_min.append(row['MutSpec_median'] - row['MutSpec_q05']) + yerr_max.append(row['MutSpec_q95'] - row['MutSpec_median']) + + yerr = [yerr_min, yerr_max] + x = list(range(len(order))) + ax.errorbar(x, ymed, yerr=yerr, fmt=".", color="gray", elinewidth=0.7, + capsize=(2 if is_192 else 4)) + else: + sns.barplot( + ms, x='Mut', y=spectra_col, hue='Mut', legend=False, width=bar_width, + order=order, palette=palette, err_kws={'linewidth': 1}, ax=ax, **kwargs) + + elif style == "box": + sns.boxplot( + ms, x='Mut', y=spectra_col, hue='Mut', legend=False, + order=order, palette=palette, ax=ax, **kwargs, + ) + else: + raise NotImplementedError + + ax.grid(axis="y", alpha=.7, linewidth=0.5) + ax.set_title(title, fontsize=titlesize, fontname=fontname) + ax.set_ylabel(ylabel if ylabel else "", fontsize=ylabelsize, fontname=fontname) + ax.set_xlabel("") + ax.set_xlim(-0.5, len(order) - 0.5) + + if is_192: + order_styled = ["" if x[0] == "_" else x for x in order] + ax.set_xticks(range(len(order_styled))) + ax.set_xticklabels(order_styled, rotation=tick_rotation, fontsize=ticksize, fontname=fontname) + else: + plt.xticks(fontsize=ticksize, fontname=fontname) + + if savepath is not None: + plt.savefig(savepath, dpi=dpi, bbox_inches="tight") + if show: + plt.show() + else: + plt.close() + return ax + + +def plot_mutspec12( + mutspec: pd.DataFrame, + spectra_col="MutSpec", + title="Spectrum", + ylabel=None, + figsize=(6, 4), + style="bar", + savepath=None, + fontname=None, + ticksize=8, + titlesize=14, + ylabelsize=12, + show=True, + dpi=300, + ax=None, + **kwargs, + ): + """ + Plot a 12-component mutational spectrum. + + A convenience wrapper around :func:`plot_mutspec` for 12-component + (SBS12) spectra. + + Arguments + --------- + mutspec: pd.DataFrame + Table containing at least a ``'Mut'`` column with 12-component SBS + codes and a column named *spectra_col* with spectrum values. + spectra_col: str + Column name used for the y-axis values. + title: str + Plot title. + ylabel: str or None + Y-axis label. Defaults to an empty string when ``None``. + figsize: tuple + Figure size ``(width, height)`` in inches. + style: str + ``'bar'`` for a bar plot or ``'box'`` for a box plot. + savepath: str or None + File path to save the figure. If ``None`` the figure is not saved. + fontname: str or None + Font family for all text elements. + ticksize: int + Font size for tick labels. + titlesize: int + Font size for the title. + ylabelsize: int + Font size for the y-axis label. + show: bool + If ``True`` call ``plt.show()``; otherwise close the figure. + dpi: int + Resolution for saved figures. + ax: matplotlib.axes.Axes or None + Axes to draw on. A new figure is created when ``None``. + **kwargs + Additional keyword arguments forwarded to the underlying seaborn + plot function. + + Return + ------ + ax: matplotlib.axes.Axes + The axes containing the plot. + """ + return plot_mutspec( + mutspec=mutspec, + spectra_col=spectra_col, + title=title, + ylabel=ylabel, + figsize=figsize, + style=style, + sbs_kind=12, + sbs_order=ordered_sbs12, + savepath=savepath, + fontname=fontname, + ticksize=ticksize, + titlesize=titlesize, + ylabelsize=ylabelsize, + show=show, + dpi=dpi, + ax=ax, + **kwargs, + ) + + +def plot_mutspec192( + mutspec192: pd.DataFrame, + spectra_col="MutSpec", + title="Mutational spectrum", + ylabel=None, + figsize=(24, 8), + style="bar", + sbs_order=ordered_sbs192_kp, + labels_style="cosmic", + savepath=None, + fontname=None, + ticksize=6, + titlesize=16, + ylabelsize=16, + bar_width=0.6, + show=True, + dpi=300, + ax=None, + **kwargs, + ): + """ + Plot a barplot of a 192-component mutational spectrum. + + Arguments + --------- + mutspec192: pd.DataFrame + Table containing 192-component mutational spectrum for one or many + species; all substitution types must be present in the table. + spectra_col: str + Column name used for the y-axis values. + title: str + Title on the plot. + ylabel: str or None + Y-axis label. + figsize: tuple + Figure size ``(width, height)`` in inches. + style: str + ``'bar'`` for a bar plot or ``'box'`` for a box plot. + sbs_order: list or None + Custom ordering of SBS192 labels on the x-axis. Defaults to + COSMIC ordering. + labels_style: str + Style for tick labels. One of ``'cosmic'``/``'long'`` (full COSMIC + strings) or ``'kk'`` (compact KK-style labels). + savepath: str or None + Path to output plot file. If ``None`` no image is written. + fontname: str or None + Font family for all text elements. + ticksize: int + Font size for tick labels. + titlesize: int + Font size for the title. + ylabelsize: int + Font size for the y-axis label. + bar_width: float + Width of bars. + show: bool + If ``True`` call ``plt.show()``; otherwise close the figure. + dpi: int + Resolution for saved figures. + ax: matplotlib.axes.Axes or None + Axes to draw on. A new figure is created when ``None``. + **kwargs + Additional keyword arguments forwarded to the underlying seaborn + plot function. + + Return + ------ + ax: matplotlib.axes.Axes + The axes containing the plot. + """ + # delegate to generic plotter for 192-component spectra + return plot_mutspec( + mutspec=mutspec192, + spectra_col=spectra_col, + title=title, + ylabel=ylabel, + figsize=figsize, + style=style, + sbs_kind=192, + sbs_order=sbs_order, + labels_style=labels_style, + savepath=savepath, + fontname=fontname, + ticksize=ticksize, + titlesize=titlesize, + ylabelsize=ylabelsize, + bar_width=bar_width, + show=show, + dpi=dpi, + ax=ax, + **kwargs, + ) + + +# def plot_mutspec192kk(mutspec192: pd.DataFrame, ylabel="MutSpec", title="Mutational spectrum", show=True, figsize=(24, 6), filepath=None): + # from .sbs_orders import ordered_sbs192_kk + # ms192 = mutspec192.copy() +# ms192["long_lbl"] = ms192.Mut.str.get(2) + ms192.Mut.str.get(4) + ": " + ms192.Mut.str.get(0) + ms192.Mut.str.get(2) + ms192.Mut.str.get(-1) +# fig = plt.figure(figsize=figsize) +# ax = fig.add_subplot(111) +# ax.grid(axis="y", alpha=.7, linewidth=0.5) +# order = _prepare_nice_labels(ordered_sbs192_kk, True) +# sns.barplot( +# x="long_lbl", y=ylabel, data=ms192, +# order=order, +# errwidth=1, ax=fig.gca(), +# ) +# plt.xticks(rotation=90, fontsize=7, fontname=None) +# ax.set_title(title) +# ax.set_xlabel("") +# ax.set_ylabel("Mutational spectrum") + +# def _coloring192kk(): +# colors = "red yellow lime blue".split() +# while True: +# for clr in colors: +# yield clr + +# # map colors to bars +# clrs_iterator = _coloring192kk() +# for bar, sbs in zip(ax.patches, order): +# if len(sbs): +# bar.set_color(next(clrs_iterator)) +# bar.set_alpha(alpha=0.9) +# bar.set_width(0.3) +# if filepath is not None: +# plt.savefig(filepath, dpi=300, bbox_inches="tight") +# if show: +# plt.show() +# else: +# plt.close() diff --git a/pymutspec/io/__init__.py b/src/pymutspec/io/__init__.py similarity index 100% rename from pymutspec/io/__init__.py rename to src/pymutspec/io/__init__.py diff --git a/pymutspec/io/auxiliary.py b/src/pymutspec/io/auxiliary.py similarity index 56% rename from pymutspec/io/auxiliary.py rename to src/pymutspec/io/auxiliary.py index 2818112..10d1ab3 100644 --- a/pymutspec/io/auxiliary.py +++ b/src/pymutspec/io/auxiliary.py @@ -5,18 +5,37 @@ def load_scheme(path: str) -> Dict[str, str]: """ + TODO deprecate parse files like scheme_birds_genes.nex (just separated genes) return dict(charset_lbl: gene_fp) """ with open(path) as handle: raw_file = handle.read() - charsets = re.findall("charset\s(\w+)\s?=\s?([\w_\.]+)(\s?:.+)?;", raw_file) + charsets = re.findall(r"charset\s(\w+)\s?=\s?([\w_\.]+)(\s?:.+)?;", raw_file) scheme = {i: os.path.basename(fp) for i, (_, fp, _) in enumerate(charsets, 1)} return scheme def get_aln_files(path: str): + """ + Return the set of ``*.fna`` alignment files in a directory. + + Arguments + --------- + path: str + Path to a directory containing alignment files. + + Return + ------ + files: set of str + Absolute paths to all ``*.fna`` files found in *path*. + + Raises + ------ + AssertionError + If *path* is not a directory. + """ assert os.path.isdir(path), "path is not directory" raw_files = os.listdir(path) files = set( diff --git a/pymutspec/io/gb.py b/src/pymutspec/io/gb.py similarity index 100% rename from pymutspec/io/gb.py rename to src/pymutspec/io/gb.py diff --git a/pymutspec/io/states.py b/src/pymutspec/io/states.py similarity index 99% rename from pymutspec/io/states.py rename to src/pymutspec/io/states.py index 9eb7417..66da220 100644 --- a/pymutspec/io/states.py +++ b/src/pymutspec/io/states.py @@ -5,7 +5,6 @@ from collections import defaultdict from typing import List -import tqdm import numpy as np import pandas as pd from Bio import SeqIO @@ -269,7 +268,7 @@ def _prepare_db(self, path_states, rewrite=False): con.close() raise ValueError(f"Inappropriate type of table, expected another columns order,\ngot {repr(header)}") - for line in tqdm.tqdm(handle, total=8652300): # TODO estimate total + for line in handle: row = line.strip().split() query = "INSERT INTO states VALUES ('{}',{},{},'{}',{},{},{},{})".format(*row) cur.execute(query) diff --git a/pymutspec/utils/__init__.py b/src/pymutspec/utils/__init__.py similarity index 100% rename from pymutspec/utils/__init__.py rename to src/pymutspec/utils/__init__.py diff --git a/pymutspec/utils/configs/log_settings.yaml b/src/pymutspec/utils/configs/log_settings.yaml similarity index 100% rename from pymutspec/utils/configs/log_settings.yaml rename to src/pymutspec/utils/configs/log_settings.yaml diff --git a/pymutspec/utils/custom_profile.py b/src/pymutspec/utils/custom_profile.py similarity index 100% rename from pymutspec/utils/custom_profile.py rename to src/pymutspec/utils/custom_profile.py diff --git a/pymutspec/utils/logging.py b/src/pymutspec/utils/logging.py similarity index 86% rename from pymutspec/utils/logging.py rename to src/pymutspec/utils/logging.py index 1ca8118..87a19b3 100644 --- a/pymutspec/utils/logging.py +++ b/src/pymutspec/utils/logging.py @@ -6,6 +6,7 @@ DEFAULT_PATH_TO_LOGCONF = os.path.join(os.path.dirname(__file__), "configs/log_settings.yaml") +# TODO remove config file requirement and write directly here the logging configuration, as it is not very complex and does not need to be changed by the user. def load_logger(path=None, stream_level: str = None, filename=None): path = path or DEFAULT_PATH_TO_LOGCONF diff --git a/tests/conftest.py b/tests/conftest.py index 5e41387..d1d7d03 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,5 +1,5 @@ import pytest -from ete3 import PhyloTree +from pymutspec.annotation.phylo_tree import Tree from pymutspec.annotation import CodonAnnotation from pymutspec.io import GenesStates @@ -30,6 +30,6 @@ def coda(): @pytest.fixture -def tree_rooted() -> PhyloTree: - t = PhyloTree(path_to_tree_rooted, format=1) +def tree_rooted() -> Tree: + t = Tree(path_to_tree_rooted, format=1) return t diff --git a/tests/data/hum_cytb_ms12syn.csv b/tests/data/hum_cytb_ms12syn.csv new file mode 100644 index 0000000..ff33cbe --- /dev/null +++ b/tests/data/hum_cytb_ms12syn.csv @@ -0,0 +1,13 @@ +Mut,ObsNum,ExpNum,MutSpec,MutSpec_q05,MutSpec_median,MutSpec_q95 +T>G,181.7552321,766,0.003321943837644736,0.0030302685019133195,0.0033655975274012708,0.0037144405770718858 +T>C,8568.3304317,1308,0.09171124651647472,0.09017670431789103,0.09199465631933745,0.09401043196216709 +T>A,218.7251507,766,0.003997643743790643,0.003557553220779003,0.003992496970327321,0.004395078358680384 +G>T,470.097816,817,0.008055644844403928,0.007449770131358198,0.008006817828481769,0.008634458815654725 +G>C,102.863887,817,0.0017626862171742108,0.001500023954426661,0.0017984243088434693,0.0019944737206306186 +G>A,7705.4716993,1872,0.05762720008836421,0.056465194076136406,0.05774179760095588,0.058930181109313105 +C>T,7283.0216979,178,0.5728299747121062,0.5657500478681502,0.5726199490622286,0.5786860565161442 +C>G,81.70515,96,0.011915501687026786,0.009976497611635757,0.01197035927266933,0.014285159452249107 +C>A,38.334658,96,0.005590549458272763,0.004030407508311741,0.005521033201526958,0.0069751741447752254 +A>T,74.108498,261,0.0039752242933777954,0.003345490937707674,0.003969197767434041,0.004828833299633618 +A>G,10441.8017335,623,0.2346505248216768,0.2307873034100324,0.23520798649935495,0.2386179513583951 +A>C,85.044906,261,0.0045618597796872225,0.0038520968077063837,0.0046462753844410485,0.005340070357253974 diff --git a/tests/data/hum_cytb_ms192syn.csv b/tests/data/hum_cytb_ms192syn.csv new file mode 100644 index 0000000..831d1df --- /dev/null +++ b/tests/data/hum_cytb_ms192syn.csv @@ -0,0 +1,193 @@ +Mut,ObsNum,ExpNum,MutSpec,MutSpec_q05,MutSpec_median,MutSpec_q95 +T[T>C]T,225.837289,78.0,0.0022171157814754824,0.001975810307541699,0.0022191242267137253,0.002492842244146696 +G[T>C]T,743.8630558,63.0,0.009041484967828127,0.008532840060708676,0.009170617732153552,0.00961827421863891 +C[T>C]T,225.1001361,49.0,0.0035177664430140585,0.003210989854349574,0.0034783259286398933,0.003996531420154894 +A[T>C]T,207.4389675,34.0,0.004671956056706138,0.004332664773455997,0.004747243140433576,0.005195516947201229 +T[G>A]T,634.0125501,107.0,0.004537340808465311,0.004261272038020317,0.004501346583365177,0.004887694720700908 +G[G>A]T,503.773226,114.0,0.0033839002709793262,0.0031198853056285366,0.003379681808229254,0.0035865176125100995 +C[G>A]T,286.021821,52.0,0.004211949426023189,0.0038552789848513725,0.0042208058794178886,0.004527434054801695 +A[G>A]T,579.3560587,173.0,0.002564405939326343,0.002418371059349158,0.0025726976646365082,0.0027243659060618326 +T[C>T]T,165.737556,4.0,0.03172840665971748,0.027683838496496475,0.031645462349245236,0.03578971565313821 +G[C>T]T,559.3476445,13.0,0.032947751775892625,0.030621457027877188,0.03289167772309278,0.03493799810188874 +C[C>T]T,272.8483878,5.0,0.04178676125556816,0.03819866233358979,0.041822868107780256,0.04551769689071283 +A[C>T]T,362.0823976,10.0,0.02772647994249007,0.024992897081607356,0.027599459754102343,0.029753233447372653 +T[A>G]T,693.6470222,31.0,0.017134212807302898,0.016301572334896562,0.017154353315224302,0.01827838089143394 +G[A>G]T,356.856544,17.0,0.016074299946715852,0.015019685909069401,0.016119187482968835,0.017405724744954033 +C[A>G]T,352.825472,27.0,0.010006529757222752,0.009124718528425143,0.010124277972152925,0.010974261587950054 +A[A>G]T,938.2515079,45.0,0.0159659248755608,0.015293594767596619,0.01606237821819738,0.016842160375266247 +T[T>G]G,17.608576,100.0,0.00013483776966677157,8.600078164944679e-05,0.00013905210553017636,0.0001918286216688674 +G[T>G]G,45.272343,70.0,0.0004952474584550715,0.0003704700509081391,0.0005171779189964996,0.0006266072274809163 +C[T>G]G,25.2498211,92.0,0.00021016374669347067,0.00014314985110963068,0.00020898686897961692,0.0002650171291308432 +A[T>G]G,26.947398,78.0,0.0002645510918061935,0.00016816549798570786,0.00025946818206696996,0.0003437900535685752 +T[T>C]G,208.305724,100.0,0.0015951022520493477,0.0014417083431685338,0.0016075772160862485,0.0018253132889487589 +G[T>C]G,414.461852,70.0,0.00453391994378515,0.004231190907872306,0.004554176832858656,0.004891315446761697 +C[T>C]G,373.235741,92.0,0.003106581286964993,0.0028307875335934498,0.003108776718296051,0.0033703104933891477 +A[T>C]G,284.071224,78.0,0.0027888166590303723,0.0025519987703585224,0.0028232430418318874,0.003067485970410502 +T[T>A]G,21.157447,100.0,0.0001620132692912208,0.00010379022263246286,0.00016329080554376008,0.00022347143669481571 +G[T>A]G,37.271772,70.0,0.00040772686218420144,0.0003031881805221618,0.0004089445394889305,0.0005151800017809984 +C[T>A]G,44.0900787,92.0,0.0003669782884759523,0.0002982251211602397,0.00037232334023634407,0.00045928232612988204 +A[T>A]G,9.697993,78.0,9.520825114464935e-05,5.461480651832759e-05,9.492436923344008e-05,0.00015074608357267224 +T[G>T]G,60.723682,156.0,0.0002980717539332016,0.0002300316349241147,0.00029736043117258064,0.00034880465395922355 +G[G>T]G,57.504451,172.0,0.00025601200968818067,0.00020556167567748328,0.0002547463628818882,0.0003118233445834013 +C[G>T]G,57.276232,66.0,0.0006645349521041727,0.0005172115214467916,0.0006684329166075174,0.000792875065786265 +A[G>T]G,47.867325,99.0,0.00037024677794481855,0.0002873239827381424,0.00038017886953217977,0.0004654386550925347 +T[G>C]G,22.691622,156.0,0.00011138539934270163,7.713595312793233e-05,0.0001081820130145039,0.00014533680867210235 +G[G>C]G,5.792486,172.0,2.5788368659508653e-05,8.86248629362575e-06,2.6165287314265367e-05,4.303148588437853e-05 +C[G>C]G,4.863435,66.0,5.642694066863122e-05,2.1210929202246175e-05,6.654306088442564e-05,9.072510688527154e-05 +A[G>C]G,1.934531,99.0,1.4963315154635606e-05,0.0,1.509926082280431e-05,3.023916407093204e-05 +T[G>A]G,813.3976613,156.0,0.003992690488512459,0.0037369877877706592,0.004003031026222983,0.0042342567225479985 +G[G>A]G,598.9462925,172.0,0.002666531744442007,0.0025114215558014845,0.002672355894584124,0.002844737175322337 +C[G>A]G,263.1038041,66.0,0.0030526043308159506,0.0027489461617277505,0.0030460958623152587,0.0033817775025172313 +A[G>A]G,829.415209,269.0,0.0023610600831221536,0.0022440719484500787,0.002394092993284016,0.0025066932102064694 +T[C>T]G,431.117285,6.0,0.05502138387396565,0.05149809554528611,0.054910448962683954,0.05880435012032357 +G[C>T]G,420.8360214,11.0,0.029295947260437064,0.027648012755758393,0.029490968611025146,0.03153671461998816 +C[C>T]G,243.80892,5.0,0.0373393634984048,0.03364454221201829,0.037433625273193404,0.04041544157253797 +A[C>T]G,421.5830887,5.0,0.0645654974139155,0.0600409185047244,0.06360259422988931,0.06867317521374819 +T[C>G]G,8.197484,6.0,0.0010462046632268326,0.0005076812508326178,0.0010032771927452927,0.0016223778184755872 +G[C>G]G,6.870761,11.0,0.00047829900878125687,0.00026006481523908914,0.0005434393751966939,0.0008212575200010919 +C[C>G]G,2.985567,5.0,0.0004572399215822043,0.0,0.0005985580273187737,0.0009156347461015808 +A[C>G]G,10.890438,5.0,0.0016678718036191643,0.000609156779247765,0.0015316873458950198,0.0024382505366824627 +T[C>A]G,9.912821000000001,6.0,0.0012651247084999344,0.0007408195668701431,0.0012703308676985658,0.0017854560480248776 +G[C>A]G,8.928319,11.0,0.0006215332082985949,0.0004106752461818704,0.0005578016708168716,0.0009609871728738014 +C[C>A]G,4.944424,5.0,0.0007572390912778606,0.0002975935870570207,0.0007537911738529204,0.0012239860473431334 +A[C>A]G,0.978598,5.0,0.00014987239367949268,0.0,0.0002962893368695542,0.00030248154584007714 +T[A>T]G,14.964729,67.0,0.00017103357672816003,9.036068413156074e-05,0.00018006941706747684,0.0002480558060012857 +G[A>T]G,9.492566,38.0,0.00019128783548942995,0.00011206436377759296,0.0001946686164689255,0.00027144585547935493 +C[A>T]G,11.795249,25.0,0.00036128874199122473,0.00023447726661354583,0.00036127838167290897,0.000537477069221037 +A[A>T]G,3.570109,23.0,0.00011886143360343728,0.0,0.00012827208533208723,0.0002374886368914705 +T[A>G]G,1131.7972819,67.0,0.012935438874607526,0.012321559409714994,0.012973896410591882,0.013440432601207522 +G[A>G]G,982.2593004,38.0,0.019793831874635136,0.01893951901070376,0.019854240147457264,0.020929166895133137 +C[A>G]G,284.6274745,25.0,0.008718145941492578,0.008050091735998728,0.008709084009323582,0.009428723732586894 +A[A>G]G,1284.9548426,65.0,0.015137767581040013,0.014521849975296843,0.015131510054538137,0.015922943012175297 +T[A>C]G,28.632415,67.0,0.00032724310261916676,0.00022213376951617288,0.0003371884660354504,0.0004243594093422037 +G[A>C]G,10.709541,38.0,0.00021581150101830264,0.00011547053466484642,0.0001983986219366196,0.00028035211351888384 +C[A>C]G,13.71839,25.0,0.0004201945940475692,0.00029052244126025276,0.0004216670250514258,0.0006041774260797472 +A[A>C]G,2.972724,23.0,9.897239449757542e-05,0.0,0.00012916473066924405,0.00019895726155374363 +T[T>G]C,0.986019,26.0,2.9040177050673066e-05,0.0,5.7455468431388925e-05,5.866671132463835e-05 +G[T>G]C,4.848855,20.0,0.00018565067204684064,7.337980677178279e-05,0.0001514750546993472,0.00029868829819521166 +C[T>G]C,3.884451,27.0,0.00011016742725982572,0.0,0.0001076763140496408,0.00021637142767551068 +A[T>G]C,3.978074,16.0,0.00019038827505464266,8.892044195665069e-05,0.00019015905398598592,0.00037833747079648987 +T[T>C]C,253.362957,44.0,0.004409382424850326,0.0038906933264980263,0.004477284309844817,0.0048966850443282915 +G[T>C]C,643.4151096,43.0,0.0114580344587543,0.010815486224311463,0.011443283964228815,0.012241913352233909 +C[T>C]C,722.598295,53.0,0.01044018963096977,0.00988524977962279,0.0104831757768435,0.011072296414689984 +A[T>C]C,500.235687,38.0,0.010080414695120157,0.00926448923779747,0.010037047347291023,0.010791698947931748 +T[T>A]C,5.890783,26.0,0.0001734950150931118,5.7671088553484536e-05,0.00017294048177653038,0.000290757634131048 +G[T>A]C,14.657394,20.0,0.0005611953845918943,0.00036818733490811905,0.0005950246224806859,0.0008276893865533306 +C[T>A]C,7.864508,27.0,0.00022304634889829154,0.00010990860177502958,0.0002220666791251434,0.0003385412840837009 +A[T>A]C,2.931861,16.0,0.00014031713801452152,0.0,9.629478646755284e-05,0.000281389583756981 +T[G>T]C,7.843369,30.0,0.00020020214063751218,9.960703988844389e-05,0.00020115167350335416,0.0003034977158672872 +G[G>T]C,10.874488,32.0,0.00026022329051978555,9.665582992909289e-05,0.00023914254973737245,0.0003791378397266956 +C[G>T]C,8.804457,17.0,0.0003965892879519731,0.00018114340643531751,0.0004345623275881145,0.0006147276890846324 +A[G>T]C,22.435219,34.0,0.000505287692827427,0.0003175288375645881,0.000503759935699351,0.000700148205831276 +T[G>C]C,6.993586,30.0,0.0001785114136454037,5.115373214176811e-05,0.00015415466014141925,0.0002567825234980084 +G[G>C]C,9.747783,32.0,0.00023326157218002603,9.592475482748365e-05,0.0002341388453151181,0.00033376076578170053 +C[G>C]C,17.653522,17.0,0.0007951879053784342,0.0005183774726578096,0.0007867069913877687,0.0010994553666122807 +A[G>C]C,9.884877,34.0,0.00022262794462638848,9.032725278668566e-05,0.00022135789314145385,0.00031315414779034496 +T[G>A]C,158.597952,40.0,0.0030361617715989205,0.002650894366397395,0.003043146749448885,0.00339632611177121 +G[G>A]C,204.4182929,57.0,0.0027462003974678927,0.002480799654665761,0.0027586606655126454,0.0031134396633379436 +C[G>A]C,120.940841,22.0,0.0042095689525128,0.003574294944761012,0.004125259835590693,0.004686630860229022 +A[G>A]C,170.4765061,60.0,0.002175708005681226,0.001926168645872658,0.00215783160263713,0.002408268221162769 +T[C>T]C,118.575253,3.0,0.03026635553833383,0.027179220940995515,0.030600728455620767,0.034389555877832295 +G[C>T]C,706.999006,9.0,0.06015387623600408,0.05720007419181876,0.060004889500267604,0.064176742621949 +C[C>T]C,428.754575,8.0,0.04103988209170388,0.03821445474987017,0.04092095772553621,0.044218388592254014 +A[C>T]C,430.121043,8.0,0.04117067879656024,0.03819816957374852,0.04112998791907074,0.04420943901326682 +T[C>G]C,2.978534,3.0,0.0007602713614030038,0.0,0.0009976663463531056,0.001518632724986161 +G[C>G]C,4.935353,5.0,0.0007558498666084187,0.0003038830174752183,0.0009028108748568836,0.0012402565760439893 +C[C>G]C,1.052533,5.0,0.00016119554723865924,0.0,0.00027601820685951584,0.0003245327453995804 +A[C>G]C,5.889812,6.0,0.0007516878081042132,0.0002495969477318147,0.0007537162967471592,0.0012613588217090575 +T[C>A]C,2.802593,3.0,0.0007153623881978613,0.0,0.0009291324962208975,0.001432900648708596 +C[C>A]C,2.9292540000000002,5.0,0.0004486155793034818,0.0,0.00030854965260457734,0.0008960595269220498 +A[C>A]C,0.992824,6.0,0.0001267092559818985,0.0,0.0,0.0002563760655710445 +G[A>T]C,6.870869,7.0,0.0007516245425220733,0.0002151905911336587,0.0006556311957602096,0.0012688015184293912 +C[A>T]C,4.119918,4.0,0.000788707380894949,0.0,0.0008156178092898051,0.0012043664094491086 +T[A>G]C,159.6905,6.0,0.020380515018147582,0.017705846824007467,0.020281836244572957,0.022663193261630736 +G[A>G]C,467.2675804,13.0,0.027523877866157974,0.025235006105909393,0.027340898345502955,0.02911470334188694 +C[A>G]C,161.675864,7.0,0.017686197672501243,0.015558175336945205,0.017713028909460395,0.01960155433449545 +A[A>G]C,266.404995,10.0,0.02039997746205453,0.018262045379882737,0.020520304174467076,0.022250686072974087 +T[A>C]C,0.956961,5.0,0.00014655868469782383,0.0,0.0002905412336779,0.0002969511568177134 +G[A>C]C,2.821708,7.0,0.00030867492665496526,0.0,0.00021774240410470062,0.0006183768185273502 +C[A>C]C,6.869567,4.0,0.0013150936976057221,0.0007348921190191751,0.001166261398785708,0.0022424698026251217 +T[T>G]A,7.823672,99.0,6.051496192229448e-05,2.9921597046662096e-05,6.086008812551125e-05,0.00010524172555928278 +G[T>G]A,13.750854,70.0,0.00015042463110616464,8.690004110149684e-05,0.00015028487340815237,0.00021499360068608106 +C[T>G]A,23.492371,103.0,0.0001746533620793334,0.00011637856342334661,0.00017388078261839658,0.0002242968141845635 +A[T>G]A,7.912798,65.0,9.321891561367955e-05,4.3947220847661e-05,9.409291475382189e-05,0.00014101449266604342 +T[T>C]A,775.339811,186.0,0.0031920263223392934,0.0029988983777540024,0.003183694678902912,0.003371810990648221 +G[T>C]A,1256.2225323,121.0,0.007950025702562051,0.00751298166586639,0.007918742516572954,0.008374834179975811 +C[T>C]A,1130.6542638,158.0,0.005479741378774306,0.005221163135053452,0.005479630341556515,0.005815124712476301 +A[T>C]A,604.1877866,101.0,0.004580763733627922,0.004392647681191924,0.0045876632072664346,0.004867898724938908 +T[T>A]A,22.38784,99.0,0.0001731666773763549,0.0001064599345875037,0.00016710257214498008,0.00022787187547650654 +G[T>A]A,26.398721,70.0,0.0002887833634259778,0.00021227333615545486,0.0002976680711545821,0.0003680717330928182 +C[T>A]A,20.498513,103.0,0.00015239561017816903,0.00010248926273302304,0.0001470030926112476,0.0002055495750043466 +A[T>A]A,5.87824,65.0,6.925023974034918e-05,2.304258968051419e-05,6.923842739644108e-05,0.00011573810534774337 +T[G>T]A,32.149807,72.0,0.00034192684411097637,0.00025105673434443687,0.00035031878174974053,0.0004443319885932778 +G[G>T]A,63.41679,76.0,0.0006389663497092211,0.0005078927293686359,0.0006299628510449402,0.0007724802444649009 +C[G>T]A,34.113609000000004,16.0,0.0016326572038123308,0.0013102562927251632,0.0016811808336978282,0.0020709964367368436 +A[G>T]A,67.088387,47.0,0.0010930419158257103,0.0008942392391535157,0.001083569999137731,0.0013168595069345793 +T[G>C]A,7.831548,72.0,8.329183724628978e-05,4.1389975959063436e-05,8.303723309724842e-05,0.0001250861073622609 +G[G>C]A,10.491022,76.0,0.00010570402620598004,5.7933487370116306e-05,9.858231615351965e-05,0.00015417217499806912 +C[G>C]A,2.99209,16.0,0.00014319966242665313,0.0,0.00018917581944606746,0.000285085637282799 +A[G>C]A,1.9873850000000002,47.0,3.2379599585294543e-05,0.0,3.2431729729278775e-05,6.514140262930563e-05 +T[G>A]A,968.1455968,183.0,0.004051136802041584,0.0038327661602749253,0.00406634333913086,0.004237234023250182 +G[G>A]A,647.9080144,158.0,0.0031401007981120554,0.0029299809752859228,0.0031646053713354225,0.0033881403602611776 +C[G>A]A,386.456873,72.0,0.004110132884153437,0.0037826999494856373,0.00413565275709637,0.004403806160671112 +A[G>A]A,540.5010004,171.0,0.0024204031858283746,0.0022645492782263793,0.0024068091969417014,0.002600106380205296 +T[C>T]A,616.026096,22.0,0.021441924053258978,0.02010322959333267,0.021399703823365722,0.022987390070267802 +G[C>T]A,898.5819479,32.0,0.02150280098559275,0.020177259452098625,0.021556113637559826,0.02289817758180721 +C[C>T]A,689.9845105,21.0,0.025159810792098215,0.0233487968235622,0.025064565466720236,0.02648539603483703 +A[C>T]A,516.6179655,16.0,0.024725031086345776,0.023064722521626807,0.024681636818333026,0.026349147738864282 +T[C>G]A,8.85777,12.0,0.0005652368629076155,0.000250292642766229,0.0005085754658800226,0.000876597783199894 +G[C>G]A,22.450412,17.0,0.0010112597414364604,0.0006145174841411309,0.0009815364552548506,0.001325717744334686 +C[C>G]A,0.978184,12.0,6.242041230540226e-05,0.0,0.00012411371503898475,0.00012614457889657745 +A[C>G]A,5.618302,9.0,0.00047802421261748445,0.00016801757747471184,0.0004638921989697788,0.000793406297714787 +T[C>A]A,2.971695,12.0,0.00018963142634300127,0.0,0.00012897953328314837,0.00038014401354255554 +G[C>A]A,0.978452,17.0,4.407353934208368e-05,0.0,8.692838032807384e-05,8.907439752688437e-05 +C[C>A]A,1.9249,12.0,0.0001228327713872531,0.0,0.00012101535996144721,0.0002471476552751193 +A[C>A]A,0.970778,9.0,8.259708877813553e-05,0.0,0.0,0.0001669193553388206 +T[A>T]A,6.779489,34.0,0.00015268816209723296,8.742288735882989e-05,0.0001726717755675242,0.00026580894318875656 +G[A>T]A,8.772161,28.0,0.00023990311482270196,0.00010568248537156333,0.0002144955452911835,0.00036777723590924094 +C[A>T]A,6.772059,11.0,0.00047142799860280246,0.00014080404557925195,0.0005311266112889867,0.0008075276818984179 +A[A>T]A,0.971349,13.0,5.721623404415986e-05,0.0,0.0,0.00011544050418564817 +T[A>G]A,1357.4668297,105.0,0.009899818929863182,0.009564102540973136,0.009895920033045182,0.010310264931655974 +G[A>G]A,912.555524,75.0,0.009317198737195031,0.008888380273157413,0.009370211033418577,0.009863781808390827 +C[A>G]A,386.2031869,45.0,0.006571895719676005,0.006098621408785711,0.006553245845165633,0.007057372132973606 +A[A>G]A,705.317808,47.0,0.011491436336400674,0.01082931500159903,0.011514228185424594,0.012151290134980614 +T[A>C]A,2.97154,34.0,6.69252477876152e-05,0.0,8.793843044405342e-05,0.00013466215272094978 +G[A>C]A,4.716034,28.0,0.00012897520305541203,0.0,0.0001489386483430225,0.0002088901383634547 +C[A>C]A,10.676026,11.0,0.0007431975371466024,0.00041073486012362955,0.0006920659306553132,0.001084371503614978 +A[C>A]T,0.0,,0.0,0.0,0.0,0.0 +A[A>C]T,0.0,,0.0,0.0,0.0,0.0 +G[T>G]T,0.0,,0.0,0.0,0.0,0.0 +G[A>C]T,0.0,,0.0,0.0,0.0,0.0 +T[A>T]C,0.0,5.0,0.0,0.0,0.0,0.0 +G[T>A]T,0.0,,0.0,0.0,0.0,0.0 +T[C>A]T,0.0,,0.0,0.0,0.0,0.0 +G[C>A]C,0.0,5.0,0.0,0.0,0.0,0.0 +T[G>T]T,0.0,,0.0,0.0,0.0,0.0 +A[A>T]T,0.0,,0.0,0.0,0.0,0.0 +T[T>G]T,0.0,,0.0,0.0,0.0,0.0 +A[G>T]T,0.0,,0.0,0.0,0.0,0.0 +C[T>G]T,0.0,,0.0,0.0,0.0,0.0 +C[C>A]T,0.0,,0.0,0.0,0.0,0.0 +T[G>C]T,0.0,,0.0,0.0,0.0,0.0 +T[C>G]T,0.0,,0.0,0.0,0.0,0.0 +A[T>A]T,0.0,,0.0,0.0,0.0,0.0 +A[G>C]T,0.0,,0.0,0.0,0.0,0.0 +G[G>T]T,0.0,,0.0,0.0,0.0,0.0 +G[G>C]T,0.0,,0.0,0.0,0.0,0.0 +C[T>A]T,0.0,,0.0,0.0,0.0,0.0 +C[A>C]T,0.0,,0.0,0.0,0.0,0.0 +C[A>T]T,0.0,,0.0,0.0,0.0,0.0 +T[A>T]T,0.0,,0.0,0.0,0.0,0.0 +A[A>C]A,0.0,13.0,0.0,0.0,0.0,0.0 +G[A>T]T,0.0,,0.0,0.0,0.0,0.0 +C[G>C]T,0.0,,0.0,0.0,0.0,0.0 +G[C>G]T,0.0,,0.0,0.0,0.0,0.0 +A[A>C]C,0.0,6.0,0.0,0.0,0.0,0.0 +A[C>G]T,0.0,,0.0,0.0,0.0,0.0 +T[A>C]T,0.0,,0.0,0.0,0.0,0.0 +A[A>T]C,0.0,6.0,0.0,0.0,0.0,0.0 +T[T>A]T,0.0,,0.0,0.0,0.0,0.0 +C[C>G]T,0.0,,0.0,0.0,0.0,0.0 +A[T>G]T,0.0,,0.0,0.0,0.0,0.0 +C[G>T]T,0.0,,0.0,0.0,0.0,0.0 +G[C>A]T,0.0,,0.0,0.0,0.0,0.0 diff --git a/tests/test_codon_ann.py b/tests/test_codon_ann.py index 54edcd0..16ecd80 100644 --- a/tests/test_codon_ann.py +++ b/tests/test_codon_ann.py @@ -124,8 +124,14 @@ def test_collect_exp_mut_freqs_on_real_gene_proba(coda, states): for sbs192, x in exp_sbs192_freqs['all'].items() \ if x - exp_sbs192_freqs['syn'].get(sbs192, 0) > 0} - assert exp_sbs12_freqs["nonsyn"] == expected_sbs12_nonsyn - assert exp_sbs192_freqs["nonsyn"] == expected_sbs192_nonsyn + # compare keys exactly but allow small numerical differences in values + assert set(exp_sbs12_freqs["nonsyn"].keys()) == set(expected_sbs12_nonsyn.keys()) + for k, v in expected_sbs12_nonsyn.items(): + assert exp_sbs12_freqs["nonsyn"][k] == pytest.approx(v, abs=1e-3) + + assert set(exp_sbs192_freqs["nonsyn"].keys()) == set(expected_sbs192_nonsyn.keys()) + for k, v in expected_sbs192_nonsyn.items(): + assert exp_sbs192_freqs["nonsyn"][k] == pytest.approx(v, abs=1e-3) def test_collect_exp_muts_on_real_gene(coda, states_most_probable): @@ -139,10 +145,22 @@ def test_collect_exp_muts_on_real_gene(coda, states_most_probable): exp_freqs_ff = exp192[exp192.Label == 'syn4f'].groupby('Mut').Mut.count().to_dict() exp_freqs_nonsyn = exp192[exp192.Label == 'nonsyn'].groupby('Mut').Mut.count().to_dict() - assert exp_sbs192_freqs['all'] == exp_freqs_all - assert exp_sbs192_freqs['syn'] == exp_freqs_syn - assert exp_sbs192_freqs['ff'] == exp_freqs_ff - assert exp_sbs192_freqs['nonsyn'] == exp_freqs_nonsyn + # compare keys exactly but allow small numerical differences in summed probabilities + assert set(exp_sbs192_freqs['all'].keys()) == set(exp_freqs_all.keys()) + for k, v in exp_freqs_all.items(): + assert exp_sbs192_freqs['all'][k] == pytest.approx(v, abs=1e-3) + + assert set(exp_sbs192_freqs['syn'].keys()) == set(exp_freqs_syn.keys()) + for k, v in exp_freqs_syn.items(): + assert exp_sbs192_freqs['syn'][k] == pytest.approx(v, abs=1e-3) + + assert set(exp_sbs192_freqs['ff'].keys()) == set(exp_freqs_ff.keys()) + for k, v in exp_freqs_ff.items(): + assert exp_sbs192_freqs['ff'][k] == pytest.approx(v, abs=1e-3) + + assert set(exp_sbs192_freqs['nonsyn'].keys()) == set(exp_freqs_nonsyn.keys()) + for k, v in exp_freqs_nonsyn.items(): + assert exp_sbs192_freqs['nonsyn'][k] == pytest.approx(v, abs=1e-3) def test_collect_exp_muts_on_real_gene_proba(coda, states): @@ -158,10 +176,22 @@ def test_collect_exp_muts_on_real_gene_proba(coda, states): exp_freqs_ff = exp192[exp192.Label == 'syn4f'].groupby('Mut').Proba.sum().to_dict() exp_freqs_nonsyn = exp192[exp192.Label == 'nonsyn'].groupby('Mut').Proba.sum().to_dict() - assert exp_sbs192_freqs['all'] == exp_freqs_all - assert exp_sbs192_freqs['syn'] == exp_freqs_syn - assert exp_sbs192_freqs['ff'] == exp_freqs_ff - assert exp_sbs192_freqs['nonsyn'] == exp_freqs_nonsyn + # compare keys exactly but allow small numerical differences in summed probabilities + assert set(exp_sbs192_freqs['all'].keys()) == set(exp_freqs_all.keys()) + for k, v in exp_freqs_all.items(): + assert exp_sbs192_freqs['all'][k] == pytest.approx(v, abs=1e-3) + + assert set(exp_sbs192_freqs['syn'].keys()) == set(exp_freqs_syn.keys()) + for k, v in exp_freqs_syn.items(): + assert exp_sbs192_freqs['syn'][k] == pytest.approx(v, abs=1e-3) + + assert set(exp_sbs192_freqs['ff'].keys()) == set(exp_freqs_ff.keys()) + for k, v in exp_freqs_ff.items(): + assert exp_sbs192_freqs['ff'][k] == pytest.approx(v, abs=1e-3) + + assert set(exp_sbs192_freqs['nonsyn'].keys()) == set(exp_freqs_nonsyn.keys()) + for k, v in exp_freqs_nonsyn.items(): + assert exp_sbs192_freqs['nonsyn'][k] == pytest.approx(v, abs=1e-3) def test_extract_mutations_simple(): diff --git a/tests/test_mutspec_calc.py b/tests/test_mutspec_calc.py index 17560fd..339f201 100644 --- a/tests/test_mutspec_calc.py +++ b/tests/test_mutspec_calc.py @@ -86,7 +86,7 @@ def test_ms192_calc(mut, cxt_freqs, use_proba, lbl_id): divisor = cxt_freqs[lbl].get(cxt, 0) if divisor == 0: continue - cond = cur_mut.Mut.str.fullmatch(sbs.replace("[", "\[").replace("]", "\]")) + cond = cur_mut.Mut.str.fullmatch(sbs.replace("[", r"\[").replace("]", r"\]")) if use_proba: expected = cur_mut[cond].ProbaFull.sum() / divisor else: diff --git a/tests/test_plot_spectrum.py b/tests/test_plot_spectrum.py index 12bba4c..4607028 100644 --- a/tests/test_plot_spectrum.py +++ b/tests/test_plot_spectrum.py @@ -1,7 +1,7 @@ import pandas as pd from pymutspec.draw import plot_mutspec12, plot_mutspec192 -from pymutspec.draw.sbs_orders import ordered_sbs192_kk +from pymutspec.draw.spectra import ordered_sbs192_kk show = False diff --git a/tests/test_tree_annotation.py b/tests/test_tree_annotation.py index 8a0b76d..ffdaa85 100644 --- a/tests/test_tree_annotation.py +++ b/tests/test_tree_annotation.py @@ -1,5 +1,4 @@ -import pytest -from pymutspec.annotation import get_ingroup_root, get_tree_len, calc_phylocoefs +from pymutspec.annotation import get_ingroup_root, get_tree_height, calc_phylocoefs def test_get_ingroup_root(tree_rooted): @@ -8,21 +7,21 @@ def test_get_ingroup_root(tree_rooted): assert round(ingroup.dist, 4) == 0.0739 -def test_get_tree_len(tree_rooted): +def test_get_tree_height(tree_rooted): ingroup = get_ingroup_root(tree_rooted) - l1 = get_tree_len(tree_rooted, 'geom_mean') - l2 = get_tree_len(ingroup, 'geom_mean') + h1 = get_tree_height(tree_rooted, 'geom_mean') + h2 = get_tree_height(ingroup, 'geom_mean') - assert round(l1, 4) == 0.3435 - assert round(l2, 4) == 0.2689 + assert round(h1, 4) == 0.3435 + assert round(h2, 4) == 0.2689 def test_calc_phylocoefs(tree_rooted): phylocoefs = calc_phylocoefs(tree_rooted) ingroup = get_ingroup_root(tree_rooted) - tl = get_tree_len(ingroup, 'geom_mean') + tl = get_tree_height(ingroup, 'geom_mean') node = 'Node1' n1_d = tree_rooted.search_nodes(name=node)[0].get_closest_leaf()[1] diff --git a/tests/test_tree_vs_ete3.py b/tests/test_tree_vs_ete3.py new file mode 100644 index 0000000..e6c74fb --- /dev/null +++ b/tests/test_tree_vs_ete3.py @@ -0,0 +1,125 @@ +""" +Tests comparing the custom Tree/TreeNode implementation against ete3's +PhyloTree/PhyloNode using the same newick file. +""" +import pytest +from ete3 import PhyloTree as Ete3PhyloTree + +from pymutspec.annotation.phylo_tree import Tree +from pymutspec.annotation import iter_tree_edges + +PATH = "./tests/data/treefile_rooted.nwk" + + +@pytest.fixture(scope="module") +def ete3_tree(): + return Ete3PhyloTree(PATH, format=1) + + +@pytest.fixture(scope="module") +def custom_tree(): + return Tree(PATH, format=1) + + +# --------------------------------------------------------------------------- +# Loading +# --------------------------------------------------------------------------- + +def test_root_name(ete3_tree, custom_tree): + """Root node names must be identical.""" + assert custom_tree.name == ete3_tree.name + + +def test_root_dist(ete3_tree, custom_tree): + """Root branch lengths must be identical.""" + assert round(custom_tree.dist, 8) == round(ete3_tree.dist, 8) + + +def test_leaf_count(ete3_tree, custom_tree): + """len(tree) returns the number of leaves in both implementations.""" + assert len(custom_tree) == len(ete3_tree) + + +def test_total_node_count(ete3_tree, custom_tree): + """get_cached_content() must return the same number of nodes.""" + assert len(custom_tree.get_cached_content()) == len(ete3_tree.get_cached_content()) + + +# --------------------------------------------------------------------------- +# Node naming +# --------------------------------------------------------------------------- + +def test_all_leaf_names_match(ete3_tree, custom_tree): + """The sorted set of leaf names must be identical.""" + ete3_leaves = sorted(n.name for n in ete3_tree.iter_leaves()) + custom_leaves = sorted(n.name for n in custom_tree.iter_leaves()) + assert custom_leaves == ete3_leaves + + +def test_all_node_names_match(ete3_tree, custom_tree): + """The sorted set of all node names must be identical.""" + ete3_nodes = sorted(n.name for n in ete3_tree.traverse()) + custom_nodes = sorted(n.name for n in custom_tree.traverse()) + assert custom_nodes == ete3_nodes + + +def test_search_nodes_by_name(ete3_tree, custom_tree): + """search_nodes(name=...) must return a node with the correct name.""" + node_name = "Node1" + ete3_match = ete3_tree.search_nodes(name=node_name) + custom_match = custom_tree.search_nodes(name=node_name) + assert len(custom_match) == len(ete3_match) == 1 + assert custom_match[0].name == ete3_match[0].name + + +# --------------------------------------------------------------------------- +# Branch / edge iteration +# --------------------------------------------------------------------------- + +def test_iter_tree_edges_count(ete3_tree, custom_tree): + """iter_tree_edges must yield the same number of edges for both trees.""" + ete3_edges = list(iter_tree_edges(ete3_tree)) + custom_edges = list(iter_tree_edges(custom_tree)) + assert len(custom_edges) == len(ete3_edges) + + +def test_iter_tree_edges_node_names(ete3_tree, custom_tree): + """The sorted (ref, alt) name pairs from iter_tree_edges must be identical.""" + def edge_name_pairs(tree): + return sorted((ref.name, alt.name) for ref, alt in iter_tree_edges(tree)) + + assert edge_name_pairs(custom_tree) == edge_name_pairs(ete3_tree) + + +def test_iter_descendants_count(ete3_tree, custom_tree): + """iter_descendants must yield the same number of nodes.""" + ete3_desc = list(ete3_tree.iter_descendants()) + custom_desc = list(custom_tree.iter_descendants()) + assert len(custom_desc) == len(ete3_desc) + + +def test_iter_leaves_names(ete3_tree, custom_tree): + """iter_leaves must yield leaves with the same names (sorted).""" + ete3_leaves = sorted(n.name for n in ete3_tree.iter_leaves()) + custom_leaves = sorted(n.name for n in custom_tree.iter_leaves()) + assert custom_leaves == ete3_leaves + + +# --------------------------------------------------------------------------- +# Branch lengths +# --------------------------------------------------------------------------- + +def test_all_branch_lengths_match(ete3_tree, custom_tree): + """ + For every named node, the branch length reported by the custom tree must + match ete3's value (to 8 decimal places). + """ + ete3_dists = {n.name: n.dist for n in ete3_tree.traverse() if n.name} + custom_dists = {n.name: n.dist for n in custom_tree.traverse() if n.name} + + assert set(custom_dists.keys()) == set(ete3_dists.keys()) + for name in ete3_dists: + assert round(custom_dists[name], 8) == round(ete3_dists[name], 8), ( + f"Branch length mismatch for node '{name}': " + f"custom={custom_dists[name]}, ete3={ete3_dists[name]}" + ) diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 0000000..f4e7d49 --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,259 @@ +""" +Tests for utility functions in pymutspec.annotation and pymutspec.constants. + +Covers: +- rev_comp / lbl2lbl_id / lbl_id2lbl (annotation.auxiliary) +- node_parent / iter_tree_edges (annotation.tree) +- get_iqr_bounds / filter_outlier_branches (annotation.spectra) +- complete_sbs192_columns / collapse_sbs192 (annotation.spectra) +- get_cossim / get_eucdist (annotation.spectra) +""" + +import pytest +import pandas as pd + +from pymutspec.annotation import ( + rev_comp, lbl2lbl_id, lbl_id2lbl, + node_parent, iter_tree_edges, + get_iqr_bounds, filter_outlier_branches, + complete_sbs192_columns, collapse_sbs192, + get_cossim, get_eucdist, +) +from pymutspec.constants import possible_sbs192, possible_sbs12 + + +# --------------------------------------------------------------------------- +# rev_comp +# --------------------------------------------------------------------------- + +class TestRevComp: + def test_pyrimidine_mutation_unchanged_context(self): + # C>A in context A_C = A[C>A]C + # rev-comp: swap flanks (C,A) → C,A → complement → G,T + # middle: [C>A] complement → [G>T] + assert rev_comp("A[C>A]C") == "G[G>T]T" + + def test_involution(self): + """Applying rev_comp twice returns the original string.""" + for sbs in ["A[C>A]C", "T[G>T]A", "C[A>G]T", "G[T>C]G"]: + assert rev_comp(rev_comp(sbs)) == sbs + + def test_all_possible_sbs192_round_trip(self): + """rev_comp is an involution on every element of possible_sbs192.""" + for sbs in possible_sbs192: + assert rev_comp(rev_comp(sbs)) == sbs + + +# --------------------------------------------------------------------------- +# lbl_id2lbl / lbl2lbl_id +# --------------------------------------------------------------------------- + +class TestLabelConversions: + @pytest.mark.parametrize("lbl_id,lbl", [(0, "all"), (1, "syn"), (2, "ff")]) + def test_lbl_id2lbl(self, lbl_id, lbl): + assert lbl_id2lbl(lbl_id) == lbl + + def test_lbl_id2lbl_invalid(self): + with pytest.raises(NotImplementedError): + lbl_id2lbl(99) + + @pytest.mark.parametrize("lbl,lbl_id", [ + ("all", 0), ("syn", 1), ("syn_c", 1), ("ff", 2), ("syn4f", 2), + ]) + def test_lbl2lbl_id(self, lbl, lbl_id): + assert lbl2lbl_id(lbl) == lbl_id + + def test_lbl2lbl_id_invalid(self): + with pytest.raises(NotImplementedError): + lbl2lbl_id("nonsyn") + + def test_round_trip(self): + for lbl_id in (0, 1, 2): + assert lbl2lbl_id(lbl_id2lbl(lbl_id)) == lbl_id + + +# --------------------------------------------------------------------------- +# node_parent / iter_tree_edges +# --------------------------------------------------------------------------- + +class TestTreeUtils: + def test_node_parent_root_returns_none(self, tree_rooted): + assert node_parent(tree_rooted) is None + + def test_node_parent_leaf(self, tree_rooted): + leaf = next(tree_rooted.iter_leaves()) + parent = node_parent(leaf) + assert parent is not None + assert leaf in parent.children + + def test_iter_tree_edges_count(self, tree_rooted): + """Number of edges equals number of nodes minus one (tree property).""" + n_nodes = sum(1 for _ in tree_rooted.traverse()) + edges = list(iter_tree_edges(tree_rooted)) + assert len(edges) == n_nodes - 1 + + def test_iter_tree_edges_root_not_yielded_as_alt(self, tree_rooted): + """The tree root must never appear as an alt_node.""" + root_name = tree_rooted.name + for _, alt in iter_tree_edges(tree_rooted): + assert alt.name != root_name + + def test_iter_tree_edges_ref_is_parent_of_alt(self, tree_rooted): + for ref, alt in iter_tree_edges(tree_rooted): + assert node_parent(alt) is ref + + +# --------------------------------------------------------------------------- +# get_iqr_bounds / filter_outlier_branches +# --------------------------------------------------------------------------- + +class TestIQRAndOutlierFiltering: + def test_get_iqr_bounds_simple(self): + s = pd.Series([1, 2, 3, 4, 5, 100]) + lb, ub = get_iqr_bounds(s) + assert lb < 1 + assert ub > 5 + assert 100 > ub # 100 is an outlier + + def test_get_iqr_bounds_symmetric(self): + s = pd.Series(range(1, 11)) + lb, ub = get_iqr_bounds(s) + assert lb < 1 + assert ub > 10 + + def test_filter_outlier_branches_removes_high_count(self): + # Create a simple mutation table with one branch having many mutations + normal = pd.DataFrame({ + "AltNode": ["A"] * 5 + ["B"] * 5 + ["C"] * 5, + "Mut": ["A[C>A]T"] * 15, + "ProbaMut": [1.0] * 15, + }) + outlier = pd.DataFrame({ + "AltNode": ["OUTLIER"] * 100, + "Mut": ["A[C>A]T"] * 100, + "ProbaMut": [1.0] * 100, + }) + obs_df = pd.concat([normal, outlier], ignore_index=True) + filtered = filter_outlier_branches(obs_df, use_proba=True) + assert "OUTLIER" not in filtered["AltNode"].values + assert set(filtered["AltNode"].unique()).issubset({"A", "B", "C"}) + + def test_filter_outlier_branches_no_proba(self): + # Need enough normal branches for IQR to identify the outlier + n_normal = 20 + normal = pd.DataFrame({ + "AltNode": [f"N{i}" for i in range(n_normal) for _ in range(5)], + "Mut": ["A[C>A]T"] * (n_normal * 5), + "ProbaMut": [1.0] * (n_normal * 5), + }) + outlier = pd.DataFrame({ + "AltNode": ["OUTLIER"] * 1000, + "Mut": ["A[C>A]T"] * 1000, + "ProbaMut": [1.0] * 1000, + }) + obs_df = pd.concat([normal, outlier], ignore_index=True) + filtered = filter_outlier_branches(obs_df, use_proba=False) + assert "OUTLIER" not in filtered["AltNode"].values + + +# --------------------------------------------------------------------------- +# complete_sbs192_columns / collapse_sbs192 +# --------------------------------------------------------------------------- + +class TestSbs192Utilities: + def test_complete_sbs192_columns_fills_zeros(self): + partial = pd.DataFrame( + {"A[C>A]A": [1.0], "T[G>T]T": [2.0]}, + ) + complete = complete_sbs192_columns(partial) + assert list(complete.columns) == possible_sbs192 + assert complete.shape == (1, 192) + # Original values preserved + assert complete["A[C>A]A"].iloc[0] == 1.0 + assert complete["T[G>T]T"].iloc[0] == 2.0 + # Missing columns filled with 0 + assert complete["A[A>C]A"].iloc[0] == 0.0 + + def test_complete_sbs192_columns_already_complete(self): + data = {sbs: [1.0] for sbs in possible_sbs192} + df = pd.DataFrame(data) + result = complete_sbs192_columns(df) + assert list(result.columns) == possible_sbs192 + assert result.shape == (1, 192) + + def test_collapse_sbs192_to_12_shape(self): + data = {sbs: [1.0] for sbs in possible_sbs192} + df = pd.DataFrame(data) + result = collapse_sbs192(df, to=12) + assert list(result.columns) == possible_sbs12 + assert result.shape == (1, 12) + + def test_collapse_sbs192_to_12_sum_preserved(self): + """Sum of 12-component spectrum equals sum of 192-component.""" + data = {sbs: [float(i)] for i, sbs in enumerate(possible_sbs192)} + df = pd.DataFrame(data) + result = collapse_sbs192(df, to=12) + assert abs(result.values.sum() - df.values.sum()) < 1e-5 + + def test_collapse_sbs192_invalid_to(self): + data = {sbs: [1.0] for sbs in possible_sbs192} + df = pd.DataFrame(data) + with pytest.raises(NotImplementedError): + collapse_sbs192(df, to=96) + + +# --------------------------------------------------------------------------- +# get_cossim / get_eucdist +# --------------------------------------------------------------------------- + +class TestDistanceMetrics: + def _make_dfs(self): + cols = possible_sbs192 + a = pd.DataFrame( + [[1.0] * 192, [0.5] * 192], + index=["x", "y"], + columns=cols, + ) + b = pd.DataFrame( + [[1.0] * 192, [1.0] * 192], + index=["x", "z"], + columns=cols, + ) + return a, b + + def test_cossim_identical_vectors(self): + cols = possible_sbs192 + df = pd.DataFrame([[1.0] * 192], index=["x"], columns=cols) + result = get_cossim(df, df) + assert abs(result["x"] - 1.0) < 1e-6 + + def test_cossim_only_common_index(self): + a, b = self._make_dfs() + result = get_cossim(a, b) + # Only "x" is common + assert list(result.index) == ["x"] + + def test_cossim_empty_on_no_common_index(self): + cols = possible_sbs192 + a = pd.DataFrame([[1.0] * 192], index=["x"], columns=cols) + b = pd.DataFrame([[1.0] * 192], index=["y"], columns=cols) + result = get_cossim(a, b) + assert len(result) == 0 + + def test_eucdist_same_vector_is_zero(self): + cols = possible_sbs192 + df = pd.DataFrame([[0.5] * 192], index=["x"], columns=cols) + result = get_eucdist(df, df) + assert abs(result["x"]) < 1e-6 + + def test_eucdist_only_common_index(self): + a, b = self._make_dfs() + result = get_eucdist(a, b) + assert list(result.index) == ["x"] + + def test_eucdist_empty_on_no_common_index(self): + cols = possible_sbs192 + a = pd.DataFrame([[1.0] * 192], index=["x"], columns=cols) + b = pd.DataFrame([[1.0] * 192], index=["y"], columns=cols) + result = get_eucdist(a, b) + assert len(result) == 0 diff --git a/tox.ini b/tox.ini index adf647d..f7ba626 100644 --- a/tox.ini +++ b/tox.ini @@ -1,2 +1,22 @@ +[tox] +minversion = 4.0 +envlist = py38,py39,py310,py311,py312,py313,py314 +isolated_build = true + [pep8] -max-line-length = 110 \ No newline at end of file +max-line-length = 110 + +[testenv] +description = Run test suite with pytest +extras = dev +deps = + pytest>=8 +commands = + pytest + +[testenv:lint] +description = run linters +skip_install = true +deps = + ruff +commands = ruff check . \ No newline at end of file