Remove S5F contents and mark experimental features (#110)

* add mutability file excerpts to docs

* make backward-compatible mutability loading

* format and lint

Author: willdumm

Diff for: Makefile

@@ -19,6 +19,8 @@ lint:
 	flake8 . --count --max-complexity=30 --max-line-length=127 --statistics
 
 docs:
+	wget -O HS5F_Mutability.csv https://bitbucket.org/kleinstein/shazam/raw/ba4b30fc6791e2cfd5712e9024803c53b136e664/data-raw/HS5F_Mutability.csv
+	wget -O HS5F_Substitution.csv https://bitbucket.org/kleinstein/shazam/raw/ba4b30fc6791e2cfd5712e9024803c53b136e664/data-raw/HS5F_Substitution.csv
 	rm -f docs/outfile docs/outtree # remove phylip dnapars output
 	make -C docs html
 	make -C docs html # need it again to get images from quickstart commands

Diff for: S5F/Mutability.csv

@@ -100,6 +100,17 @@ For example here is the top ranked tree ``gctree.out.inference.1.svg``:
 
 You will also see Python pickle files ``gctree.out.inference.[1,2,...].p`` containing a :obj:`gctree.CollapsedTree` object for each optimal tree, which can be loaded and manipulated via the API (e.g. plotted in various ways using :meth:`gctree.CollapsedTree.render`).
 
+All parsimony trees found by dnapars, as well as branching process parameters
+are saved in the file ``gctree.out.inference.parsimony_forest.p``, containing a :class:`gctree.CollapsedForest` object.
+This file may be manipulated using ``gctree infer``, instead of providing
+a dnapars ``outfile``.
+
+.. note::
+  Although described below, using mutability parsimony or isotype parsimony
+   as ranking criteria is experimental, and has not yet been shown in a careful
+   validation to improve tree inference. Only the default branching process
+   likelihood is recommended for tree ranking!
+
 Criteria other than branching process likelihoods can be used to break ties
 between trees. Providing arguments ``--isotype_mapfile`` and
 ``--idmapfile`` will allow trees to be ranked by isotype parsimony. Providing
@@ -109,15 +120,16 @@ lexicographically, first maximizing likelihood, then minimizing isotype
 parsimony and mutabilities, if such information is provided.
 Ranking priorities can be adjusted using the argument ``--ranking_coeffs``.
 
-All parsimony trees found by dnapars, as well as branching process parameters
-are saved in the file ``gctree.out.inference.parsimony_forest.p``, containing a :class:`gctree.CollapsedForest` object.
-This file may be manipulated using ``gctree infer``. For example, to find the optimal tree
+For example, to find the optimal tree
 according to a linear combination of likelihood, isotype parsimony,
 mutabilities, and alleles:
 
-.. command-output:: gctree infer gctree.out.inference.parsimony_forest.p --frame 1 --idmap idmap.txt --isotype_mapfile ../example/isotypemap.txt --mutability ../S5F/Mutability.csv --substitution ../S5F/Substitution.csv --ranking_coeffs 1 0.1 0 --outbase newranking --summarize_forest --tree_stats --verbose
+.. command-output:: gctree infer gctree.out.inference.parsimony_forest.p --frame 1 --idmap idmap.txt --isotype_mapfile ../example/isotypemap.txt --mutability ../HS5F_Mutability.csv --substitution ../HS5F_Substitution.csv --ranking_coeffs 1 0.1 0 --outbase newranking --summarize_forest --tree_stats --verbose
    :shell:
 
+The files ``HS5F_Mutability.csv`` and ``HS5F_Substitution.csv`` are a context
+sensitive mutation model which can be downloaded from the `Shazam Project <https://bitbucket.org/kleinstein/shazam/src/master/data-raw/`>_.
+
 By default, only the files listed above will be generated, with the optional argument ``--outbase`` specifying how the output files should be named.
 
 .. image:: newranking.inference.1.svg

Diff for: S5F/Substitution.csv

@@ -1,8 +1,6 @@
-r"""
-This module contains classes for simulation and inference for a binary
+r"""This module contains classes for simulation and inference for a binary
 branching process with mutation in which the tree is collapsed to nodes that
-count the number of clonal leaves of each type.
-"""
+count the number of clonal leaves of each type."""
 
 from __future__ import annotations
 
@@ -378,7 +376,8 @@ def ll(
         p: np.float64,
         q: np.float64,
     ) -> Tuple[np.float64, np.ndarray]:
-        r"""Log likelihood of branching process parameters :math:`(p, q)` given tree topology :math:`T` and genotype abundances :math:`A`.
+        r"""Log likelihood of branching process parameters :math:`(p, q)` given
+        tree topology :math:`T` and genotype abundances :math:`A`.
 
         .. math::
             \ell(p, q; T, A) = \log\mathbb{P}(T, A \mid p, q)
@@ -833,7 +832,8 @@ def support(
         weights: Optional[List[np.float64]] = None,
         compatibility: bool = False,
     ):
-        r"""Compute support from a list of bootstrap :class:`CollapsedTree` objects, and add to tree attibute.
+        r"""Compute support from a list of bootstrap :class:`CollapsedTree`
+        objects, and add to tree attibute.
 
         Args:
             bootstrap_trees_list: List of trees
@@ -867,10 +867,9 @@ def local_branching(
         self, tau=1, tau0=1, infinite_root_branch=True, nan_root_lbr=False
     ):
         r"""Add local branching statistics (Neher et al. 2014) as tree node
-        features to the ETE tree attribute.
-        After execution, all nodes will have new features ``LBI``
-        (local branching index) and ``LBR`` (local branching ratio, below Vs
-        above the node)
+        features to the ETE tree attribute. After execution, all nodes will
+        have new features ``LBI`` (local branching index) and ``LBR`` (local
+        branching ratio, below Vs above the node)
 
         Args:
             tau: decay timescale for exponential filter
@@ -994,7 +993,8 @@ def __init__(
         self.is_isotyped = False
 
     def simulate(self, p: np.float64, q: np.float64, n_trees: int):
-        r"""Simulate a forest of collapsed trees. Overwrites existing forest attribute.
+        r"""Simulate a forest of collapsed trees. Overwrites existing forest
+        attribute.
 
         Args:
             p: branching probability
@@ -1017,7 +1017,10 @@ def ll(
         q: np.float64,
         marginal: bool = False,
     ) -> Tuple[np.float64, np.ndarray]:
-        r"""Log likelihood of branching process parameters :math:`(p, q)` given tree topologies :math:`T_1, \dots, T_n` and corresponding genotype abundances vectors :math:`A_1, \dots, A_n` for each of :math:`n` trees in the forest.
+        r"""Log likelihood of branching process parameters :math:`(p, q)` given
+        tree topologies :math:`T_1, \dots, T_n` and corresponding genotype
+        abundances vectors :math:`A_1, \dots, A_n` for each of :math:`n` trees
+        in the forest.
 
         If ``marginal=False`` (the default), compute the joint log likelihood
 
@@ -1637,6 +1640,7 @@ def f(x):
 
 def _lltree(cm_counts, p: np.float64, q: np.float64) -> Tuple[np.float64, np.ndarray]:
     r"""Log likelihood of branching process parameters :math:`(p, q)`
+
     .. math::
         \ell(p, q; T, A) = \log\mathbb{P}(T, A \mid p, q)
 

Diff for: docs/quickstart.rst

@@ -592,7 +592,9 @@ def get_parser():
         default=None,
         help=(
             "Path to mutability model file. If --mutability and --substitution are both provided, "
-            "they will be used to rank trees after likelihood and isotype parsimony."
+            "they will be used to rank trees after likelihood and isotype parsimony. This shall be a csv file"
+            "with the first column containing fivemers, and the second column containing mutability scores."
+            "See a file excerpt in the documentation for :meth:`mutation_model.MutationModel`."
         ),
     )
     parser_infer.add_argument(
@@ -602,6 +604,9 @@ def get_parser():
         help=(
             "Path to substitution model file. If --mutability and --substitution are both provided, "
             "they will be used to rank trees after likelihood and isotype parsimony."
+            "This shall be a csv file with the first column containing fivemers, and the next four"
+            "columns containing targeting probabilities for bases A, C, G, and T, respectively."
+            "See a file excerpt in the documentation for :meth:`mutation_model.MutationModel`."
         ),
     )
     parser_infer.add_argument(

Diff for: gctree/branching_processes.py

@@ -12,7 +12,6 @@
 
      $ mkconfig sequence.phy > dnapars.cfg
      $ dnapars < dnapars.cfg
-
 """
 import os
 import argparse

Diff for: gctree/cli.py

@@ -22,6 +22,41 @@ class MutationModel:
         mutation_order: whether or not to mutate sequences using a context sensitive manner
                         where mutation order matters
         with_replacement: allow the same position to mutate multiple times on a single branch
+
+    Notes:
+        ``mutability_file`` shall be a csv file with the first column containing fivemers,
+        and the second column containing mutability scores.
+        An example can be found at https://bitbucket.org/kleinstein/shazam/src/master/data-raw/HS5F_Mutability.csv
+
+        For example:
+
+        .. code-block:: text
+
+            Fivemer,Mutability,...
+            TCGGG,0.03542,...
+            GCCGG,0.02241675,...
+            GCCGC,0.06789,...
+            .
+            .
+            .
+
+
+        ``substitution_file`` shall be a csv file with the first column containing fivemers,
+        and the next four columns containing targeting probabilities for bases A, C, G,
+        and T, respectively.
+        An example can be found at https://bitbucket.org/kleinstein/shazam/src/master/data-raw/HS5F_Substitution.csv
+
+        For example:
+
+        .. code-block:: text
+
+            Fivemer,A,C,G,T,...
+            AAAAA,0,0.33,0.33,0.34,...
+            AAAAC,0,0.5000,0.2500,0.2500,...
+            AAAAG,0,0.65,0.15,0.20,...
+            .
+            .
+            .
     """
 
     def __init__(
@@ -39,7 +74,10 @@ def __init__(
                 # eat header
                 f.readline()
                 for line in f:
-                    motif, score = line.replace('"', "").split()[:2]
+                    if line[0] == '"':
+                        motif, score = line.replace('"', "").split()[:2]
+                    else:
+                        motif, score = line.replace(",", " ").split()[:2]
                     self.context_model[motif] = float(score)
 
             # kmer k
@@ -48,7 +86,10 @@ def __init__(
                 # eat header
                 f.readline()
                 for line in f:
-                    fields = line.replace('"', "").split()
+                    if line[0] == '"':
+                        fields = line.replace('"', "").split()
+                    else:
+                        fields = line.replace(",", " ").split()
                     motif = fields[0]
                     if self.k is None:
                         self.k = len(motif)
@@ -69,8 +110,9 @@ def _disambiguate(sequence):
         return _sequence_disambiguations(sequence)
 
     def mutability(self, kmer: str) -> Tuple[np.float64, np.float64]:
-        r"""Returns the mutability of a central base of :math:`k`-mer, along with
-        nucleotide bias averages over ambiguous ``"N"`` nucleotide identities.
+        r"""Returns the mutability of a central base of :math:`k`-mer, along
+        with nucleotide bias averages over ambiguous ``"N"`` nucleotide
+        identities.
 
         Args:
             kmer: nucleotide :math:`k`-mer
@@ -213,7 +255,8 @@ def simulate(
         n: Optional[int] = None,
         verbose: bool = False,
     ) -> TreeNode:
-        r"""Simulate a neutral binary branching process with the mutation model, returning a :class:`ete3.Treenode` object.
+        r"""Simulate a neutral binary branching process with the mutation model,
+        returning a :class:`ete3.Treenode` object.
 
         Args:
             sequence: root nucleotide sequence

Diff for: gctree/mkconfig.py

@@ -1,10 +1,13 @@
 #!/bin/bash
+set -eu
 
 export QT_QPA_PLATFORM=offscreen
 export XDG_RUNTIME_DIR=/tmp/runtime-runner
 export MPLBACKEND=agg
 mkdir -p tests/smalltest_output
+wget -O HS5F_Mutability.csv https://bitbucket.org/kleinstein/shazam/raw/ba4b30fc6791e2cfd5712e9024803c53b136e664/data-raw/HS5F_Mutability.csv
+wget -O HS5F_Substitution.csv https://bitbucket.org/kleinstein/shazam/raw/ba4b30fc6791e2cfd5712e9024803c53b136e664/data-raw/HS5F_Substitution.csv
 gctree infer tests/small_outfile tests/abundances.csv --outbase tests/smalltest_output/gctree.infer  --root GL --frame 1 --verbose --idlabel
 gctree infer tests/small_outfile tests/abundances.csv --outbase tests/smalltest_output/gctree.infer --root GL --frame 1 --verbose --idlabel --idmapfile tests/idmap.txt --isotype_mapfile tests/isotypemap.txt
-gctree infer tests/small_outfile tests/abundances.csv --outbase tests/smalltest_output/gctree.infer --root GL --frame 1 --verbose --idlabel --mutability S5F/Mutability.csv --substitution S5F/Substitution.csv
-gctree infer tests/small_outfile tests/abundances.csv --outbase tests/smalltest_output/gctree.infer --root GL --frame 1 --verbose --idlabel --idmapfile tests/idmap.txt --isotype_mapfile tests/isotypemap.txt --mutability S5F/Mutability.csv --substitution S5F/Substitution.csv
+gctree infer tests/small_outfile tests/abundances.csv --outbase tests/smalltest_output/gctree.infer --root GL --frame 1 --verbose --idlabel --mutability HS5F_Mutability.csv --substitution HS5F_Substitution.csv
+gctree infer tests/small_outfile tests/abundances.csv --outbase tests/smalltest_output/gctree.infer --root GL --frame 1 --verbose --idlabel --idmapfile tests/idmap.txt --isotype_mapfile tests/isotypemap.txt --mutability HS5F_Mutability.csv --substitution HS5F_Substitution.csv