Fix spectrum, tree, IO, and plotting bugs; add calculate_mutrate

CHANGELOG.md

@@ -1,14 +1,37 @@
 # CHANGELOG
 
-<!-- ## 0.0.XX (2025-XX-XX)
+## 0.0.16 (2026-08-12)
 
 Features:
 
-- 
+- Added `calculate_mutrate` for unscaled observed/expected mutation rates
+- `calculate_mutspec` now keeps `RawMutSpec` (obs/exp) separately from scaled `MutSpec`
+- `CodonAnnotation` accepts sequences whose length is not divisible by 3 (incomplete last codon is ignored for codon-aware labels)
+- `Tree` can be loaded from a Newick string as well as a file path
+- `TreeNode` gained `parent` / `up` properties and `iter_edges()`
+- `calc_edgewise_spectra` accepts `nobs_cutoff` (the misspelled `nobs_cuttof` remains as an alias)
 
 Fixes:
 
--  -->
+- `_prepare_codontable` now raises `ValueError` for invalid genetic-code arguments
+- `collect_exp_muts_proba` no longer requires a `logger` attribute on `CodonAnnotation`
+- `GenomeStates` reads the states file instead of the gappy-sites path, and `path_to_gappy_sites` is optional
+- `GenesStates` passes `states_fmt` through to alignment readers and uses parameterized SQL
+- `GenesStates` rate-length mismatch raises a valid `RuntimeError` (previously passed an illegal `file=` keyword)
+- `MutSpecExtractor._derive_mutspec` now processes `Branch` objects (serial extraction was broken after the parallel refactor)
+- `plot_mutspec` no longer closes a figure the caller passed in via `ax=`
+- `filter_short_exp_seqs` actually drops IQR outliers (previously kept every node)
+- `scripts/calculate_mutspec.py` uses keyword arguments for `DataFrame.pivot` (required by pandas 3)
+- `collapse_mutspec` accepts `ObsNum`/`ExpNum` as well as `ObsFr`/`ExpFr`
+- `jackknife_spectra_sampling` / `calc_edgewise_spectra` no longer mutate caller-provided DataFrames
+- `filter_outlier_branches` falls back to `ProbaFull` when `ProbaMut` is absent
+- `read_genbank_ref` writes the detected gene qualifier instead of a leftover loop variable
+- `basic_logger` no longer attaches duplicate handlers
+- Spectrum scaling of an all-zero vector no longer produces NaNs
+- `complete_sbs192_columns` fills missing columns in one concat (avoids pandas fragmentation warning)
+- Tests for `calculate_mutspec` now actually assert `RawMutSpec` values
+
+**Full Changelog**: https://github.com/mitoclub/PyMutSpec/compare/0.0.15...0.0.16
 
 ## 0.0.15 (2026-07-16)
 

README.md

@@ -404,7 +404,7 @@ 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
+    nobs_cutoff=10,        # minimum observed mutations per branch
     scale=True,            # normalise each branch spectrum
 )
 # spectra_per_edge: DataFrame indexed by (RefNode, AltNode), 192 columns
@@ -428,7 +428,8 @@ cossim = get_cossim(spectra_per_edge.loc[["nodeA"]], spectra_per_edge.loc[["node
 | `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 |
+| `calculate_mutspec(obs, exp, ...)` | Compute a 12- or 192-component mutational spectrum (`RawMutSpec` = obs/exp, `MutSpec` = scaled) |
+| `calculate_mutrate(obs, exp, ...)` | Compute unscaled mutation rates (`MutRate` = observed / expected) |
 | `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` |
@@ -509,5 +510,5 @@ twine upload dist/*
 - [ ] 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
+- [x] Add feature to calculate mutation rate (MutRate)
+- [x] Rename some functions and variables for better readability

scripts/calculate_mutspec.py

@@ -28,9 +28,9 @@ def dump_expected(exp, path):
 
 def filter_short_exp_seqs(exp_freqs: pd.DataFrame):
     gene_len_proxy = exp_freqs.set_index(['Label', 'Node'])[possible_sbs12]\
-        .sum(1).unstack(0).iloc[:, 0]
+        .sum(axis=1).unstack(0).iloc[:, 0]
     lower_bound, upper_bound = get_iqr_bounds(gene_len_proxy)
-    used_nodes = gene_len_proxy.between(lower_bound, upper_bound).index.values
+    used_nodes = gene_len_proxy[gene_len_proxy.between(lower_bound, upper_bound)].index
     exp_freqs_flt = exp_freqs[exp_freqs.Node.isin(used_nodes)]
     return exp_freqs_flt
 
@@ -115,7 +115,7 @@ def main(
         if branches:
             raise ValueError("For branch specific spectra expected mutations required for every internal tree node")
 
-        exp_mean = exp_raw.pivot("Label", "Mut", "Count")
+        exp_mean = exp_raw.pivot(index="Label", columns="Mut", values="Count")
         exp_freqs = None
     else:
         raise RuntimeError("Expected another columns in the table {}".format(path_to_exp))

scripts/collect_mutations.py

@@ -172,8 +172,8 @@ def extract_mutspec_from_tree(self):
                     )
                     add_header["freqs"] = False
                     if self._save_exp_muts and node_expected_sbs is not None:
-                        self.dump_table(node_expected_sbs, self.handle["exp"], add_header["exp"])
-                        add_header["exp"] = False
+                        if self.dump_table(node_expected_sbs, self.handle["exp"], add_header["exp"]):
+                            add_header["exp"] = False
 
                 # summarize state frequencies over genome
                 for lbl in self.mut_labels:
@@ -212,8 +212,8 @@ def extract_mutspec_from_tree(self):
                         mutspec12["Label"] = lbl
                         mutspec12["Gene"]  = gene
                         # Dump gene mutspecs 
-                        self.dump_table(mutspec12, self.handle["ms12s"], add_header["ms12g"])
-                        add_header["ms12g"] = False
+                        if self.dump_table(mutspec12, self.handle["ms12s"], add_header["ms12g"]):
+                            add_header["ms12g"] = False
 
                         if gene_mut_df.Mut.nunique() >= self.mnum192:
                             if lbl == 'nonsyn':
@@ -233,8 +233,8 @@ def extract_mutspec_from_tree(self):
                             mutspec192["Label"] = lbl
                             mutspec192["Gene"] = gene
                             # Dump gene mutspecs 
-                            self.dump_table(mutspec192, self.handle["ms192g"], add_header["ms192g"])
-                            add_header["ms192g"] = False
+                            if self.dump_table(mutspec192, self.handle["ms192g"], add_header["ms192g"]):
+                                add_header["ms192g"] = False
 
             visited_nodes.add(ref_node.name)
 
@@ -252,8 +252,8 @@ def extract_mutspec_from_tree(self):
                 logger.warning(f"Observed too many mutations ({mut_num} > {aln_size} * 0.1) for branch ({ref_node.name} - {alt_node.name})")
 
             # dump mutations
-            self.dump_table(genome_mutations_df, self.handle["mut"], add_header["mut"])
-            add_header["mut"] = False
+            if self.dump_table(genome_mutations_df, self.handle["mut"], add_header["mut"]):
+                add_header["mut"] = False
 
             # calculate full genome mutational spectra for all labels
             if self.derive_spectra:
@@ -289,9 +289,10 @@ def extract_mutspec_from_tree(self):
                     mutspec192["Label"] = lbl
 
                     # Dump genome spectra
-                    self.dump_table(mutspec12,  self.handle["ms12"],  add_header["ms"])
-                    self.dump_table(mutspec192, self.handle["ms192"], add_header["ms"])
-                    add_header["ms"] = False
+                    if self.dump_table(mutspec12,  self.handle["ms12"],  add_header["ms12"]):
+                        add_header["ms12"] = False
+                    if self.dump_table(mutspec192, self.handle["ms192"], add_header["ms192"]):
+                        add_header["ms192"] = False
 
         logger.info(f"Processed {ei} tree edges")
         logger.info(f"Observed {total_mut_num:.3f} substitutions")
@@ -302,9 +303,12 @@ def extract_mutspec_from_tree(self):
 
     @staticmethod
     def dump_table(df: pd.DataFrame, handle, header=False):
+        if df is None or df.empty:
+            return False
         if header:
-            handle.write("\t".join(df.columns) + "\n")
+            handle.write("\t".join(map(str, df.columns)) + "\n")
         handle.write(df.to_csv(sep="\t", index=None, header=None, float_format='%g'))
+        return True
 
     def dump_expected_mutations(self, gene_exp_sbs12, gene_exp_sbs192, node, gene, handle, header=False):
         # TODO rewrite using self.dump_table

scripts/collect_mutations_parallel.py

@@ -380,7 +380,7 @@ def _derive_mutspec(self):
         total_mut_num = 0
 
         for edge_data in self.iter_branches():
-            edge_mutations = self.process_branch(*edge_data[1:])
+            edge_mutations = edge_data.process_branch()
             mut_num = edge_mutations['ProbaFull'].sum() if self.use_proba and \
                 'ProbaFull' in edge_mutations.columns else len(edge_mutations)
             total_mut_num += mut_num
@@ -390,7 +390,7 @@ def _derive_mutspec(self):
             add_header["mut"] = False
 
         self.close_handles()
-        logger.info(f"Processed {edge_data[1]} tree edges")
+        logger.info(f"Processed {edge_data.index} tree edges")
         logger.info(f"Observed {total_mut_num:.3f} substitutions")
         logger.info("Extraction of mutations from phylogenetic tree completed succesfully")
 
@@ -400,7 +400,11 @@ def _derive_mutspec_parallel(self):
         with mp.Pool(processes=self.num_processes) as pool:
             genome_mutations_lst = pool.map(Branch.process_branch, self.iter_branches())
 
-        genome_mutations = pd.concat(genome_mutations_lst)
+        frames = [df for df in genome_mutations_lst if df is not None and not df.empty]
+        if frames:
+            genome_mutations = pd.concat(frames, ignore_index=True)
+        else:
+            genome_mutations = pd.DataFrame()
 
         self.open_handles(self.outdir)
         self.dump_table(genome_mutations, self.handle["mut"], True)
@@ -450,9 +454,12 @@ def get_edges_data(self):
 
     @staticmethod
     def dump_table(df: pd.DataFrame, handle, header=False):
+        if df is None or df.empty:
+            return False
         if header:
-            handle.write("\t".join(df.columns) + "\n")
+            handle.write("\t".join(map(str, df.columns)) + "\n")
         handle.write(df.to_csv(sep="\t", index=None, header=None, float_format='%g'))
+        return True
 
     def turn_to_MAP(self, states: np.ndarray):
         if isinstance(states, pd.DataFrame):

src/pymutspec/__init__.py

@@ -1,8 +1,8 @@
 from .draw import plot_mutspec, plot_mutspec12, plot_mutspec192
 from .annotation import (
     CodonAnnotation, mutations_summary, 
-    calculate_mutspec, complete_sbs192_columns,
+    calculate_mutspec, calculate_mutrate, complete_sbs192_columns,
     rev_comp, transcriptor
 )
 
-__version__ = "0.0.15"
+__version__ = "0.0.16"

src/pymutspec/annotation/__init__.py

@@ -2,7 +2,7 @@
 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, 
+    calculate_mutspec, calculate_mutrate, collapse_mutspec, calc_edgewise_spectra, 
     complete_sbs192_columns, jackknife_spectra_sampling, collapse_sbs192, 
     get_cossim, get_eucdist, filter_outlier_branches, sample_spectrum,
     get_iqr_bounds,