Initial commit

Author: jmenglund

.gitignore

@@ -0,0 +1,53 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+
+# C extensions
+*.so
+
+# Distribution / packaging
+bin/
+build/
+develop-eggs/
+dist/
+eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+.tox/
+.coverage
+.cache
+nosetests.xml
+coverage.xml
+
+# Translations
+*.mo
+
+# Mr Developer
+.mr.developer.cfg
+.project
+.pydevproject
+
+# Rope
+.ropeproject
+
+# Django stuff:
+*.log
+*.pot
+
+# Sphinx documentation
+docs/_build/
+
+# Apple OSX files
+.DS_Store

.travis.yml

@@ -0,0 +1,18 @@
+language: python
+
+python:
+    - '2.7'
+    - '3.3'
+    - '3.4'
+    - '3.5'
+
+# whitelist
+branches:
+    only:
+        - master
+
+script:
+    py.test --cov=pandascharm.py --pep8
+
+after_success:
+  - codecov

CHANGELOG.rst

@@ -0,0 +1,10 @@
+Changelog
+=========
+
+
+v0.1.0
+------
+
+Initial release.
+
+Release date: 2016-08-??

LICENSE.txt

@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2016 Markus Englund
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

MANIFEST.in

@@ -0,0 +1,2 @@
+include *.rst
+include LICENSE.txt

README.rst

@@ -0,0 +1,103 @@
+predsim
+=======
+
+|License|
+
+``predsim`` is a command-line tool for simulating predictive
+datasets from `MrBayes <http://mrbayes.sourceforge.net>`_ output files. 
+Datasets can be simulated under the GTR + G + I substitution model or any of
+its nested variants available in MrBayes (JC69, HKY85 etc.).
+
+The script uses `Seq-Gen <http://tree.bio.ed.ac.uk/software/seqgen/>`_ for 
+simulating the DNA-sequences and builds on the third-party libraries 
+`DendroPy <http://dendropy.org>`_ and `pandas <http://pandas.pydata.org>`_.
+
+Source repository: `<https://github.com/jmenglund/predsim>`_
+
+--------------------------------
+
+.. contents:: Table of contents
+   :backlinks: top
+   :local:
+
+
+Requirements
+------------
+
+`SeqGen <http://tree.bio.ed.ac.uk/software/seqgen/>`_ must be installed on
+your system.
+
+
+Installation
+------------
+
+For most users, the easiest way is probably to install the latest version 
+hosted on `PyPI <https://pypi.python.org/>`_:
+
+.. code-block:: bash
+
+    $ pip install predsim
+
+The project is hosted at https://github.com/jmenglund/predsim and 
+can also be installed using git:
+
+.. code-block:: bash
+
+    $ git clone https://github.com/jmenglund/predsim.git
+    $ cd predsim
+    $ python setup.py install
+
+
+You may consider installing ``predsim`` and its required Python packages 
+within a virtual environment in order to avoid cluttering your system's 
+Python path. See for example the environment management system 
+`conda <http://conda.pydata.org>`_ or the package 
+`virtualenv <https://virtualenv.pypa.io/en/latest/>`_.
+
+
+Usage
+-----
+
+.. code-block:: bash
+    
+    $ predsim -h
+    usage: predsim [-h] [-V] [-l INT] [-g INT] [-c PATH] [-s INT] [-p PATH]
+                       pfile tfile [outfile]
+    
+    A command-line utility that reads posterior output of MrBayes and simulates
+    predictive datasets with SeqGen.
+    
+    positional arguments:
+      pfile                 path to a MrBayes p-file
+      tfile                 path to a MrBayes t-file
+      outfile               path to output file (default: <stdout>)
+    
+    optional arguments:
+      -h, --help            show this help message and exit
+      -V, --version         show program's version number and exit
+      -l INT, --length INT  sequence lenght (default: 1000)
+      -g INT, --gamma-cats INT
+                            number of gamma rate categories (default: continuous)
+      -c PATH, --commands-file PATH
+                            path to output file with used SeqGen commands
+      -s INT, --skip INT    number of records (trees) to skip at the beginning of
+                            the sample (default: 0)
+      -p PATH, --seqgen-path PATH
+                            path to a SeqGen executable (default: "seq-gen")
+
+
+License
+-------
+
+``predsim`` is distributed under the 
+`MIT license <https://opensource.org/licenses/MIT>`_.
+
+
+Author
+------
+
+Markus Englund, `orcid.org/0000-0003-1688-7112 <http://orcid.org/0000-0003-1688-7112>`_
+
+
+.. |License| image:: https://img.shields.io/badge/license-MIT-blue.svg
+   :target: https://raw.githubusercontent.com/jmenglund/predsim/master/LICENSE.txt

predsim.py

@@ -0,0 +1,265 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+Command-line tool for simulating predictive datasets from MrBayes' output.
+"""
+
+import argparse
+import math
+import sys
+
+import dendropy
+import pandas
+
+
+__authors__ = 'Markus Englund'
+__license__ = 'MIT'
+__version__ = '0.1.0'
+
+
+def _get_skiprows(cnt):
+    """
+    Return a list with rows to skip when
+    reading the p-file.
+
+    Parameters
+    ----------
+    cnt : int
+        Number of rows to be skipped.
+    """
+    skiprows = [0]
+    if cnt > 0:
+        skiprows.extend(range(2, cnt + 2))
+    return skiprows
+
+
+def _iterrecords(frame):
+    """Iterate over the rows of a DataFrame as dicts."""
+    for record in frame.to_dict('records'):
+        yield record
+
+
+def kappa_to_titv(kappa, piA, piC, piG, piT):
+    """Calculate transistion/transversion ratio from kappa."""
+    tot = piA + piC + piG + piT
+    if math.fabs(tot - 1.0) > 1.e-6:
+        piA = piA/tot
+        piC = piC/tot
+        piG = piG/tot
+        piT = piT/tot
+    titv = kappa * (piA * piG + piC * piT)/((piA + piG) * (piC + piT))
+    return titv
+
+
+def get_seqgen_params(mrbayes_params):
+    """
+    Adapt MrBayes parameter values for use with Seq-Gen.
+
+    Paramters
+    ---------
+    mrbayes_prams : dict
+        Parameter values from a single row in a MrBayes p-file.
+
+    Returns
+    -------
+    seqgen_params : dict
+    """
+    seqgen_params = {}
+    try:
+        seqgen_params['state_freqs'] = (
+            str(mrbayes_params['pi(A)']) + ',' +
+            str(mrbayes_params['pi(C)']) + ',' +
+            str(mrbayes_params['pi(G)']) + ',' +
+            str(mrbayes_params['pi(T)']))
+    except KeyError as ex:
+        raise KeyError(
+            'Could not find any base frequences:\n{ex}'.format(ex=str(ex)))
+    try:
+        seqgen_params['ti_tv'] = kappa_to_titv(
+            float(mrbayes_params['kappa']),
+            float(mrbayes_params['pi(A)']),
+            float(mrbayes_params['pi(C)']),
+            float(mrbayes_params['pi(G)']),
+            float(mrbayes_params['pi(T)']))
+    except KeyError:
+        pass
+    try:
+        seqgen_params['general_rates'] = (
+            str(mrbayes_params['r(A<->C)']) + ',' +
+            str(mrbayes_params['r(A<->G)']) + ',' +
+            str(mrbayes_params['r(A<->T)']) + ',' +
+            str(mrbayes_params['r(C<->G)']) + ',' +
+            str(mrbayes_params['r(C<->T)']) + ',' +
+            str(mrbayes_params['r(G<->T)']))
+    except KeyError:
+        pass
+    try:
+        seqgen_params['gamma_shape'] = str(mrbayes_params['alpha'])
+    except KeyError:
+        pass
+    try:
+        seqgen_params['prop_invar'] = str(mrbayes_params['pinvar'])
+    except KeyError:
+        pass
+    return seqgen_params
+
+
+def simulate_matrix(
+        seqgen_path, seqgen_params, tree, seq_len=1000, gamma_cats=None):
+    """
+    Simulate a predictive datasets with Seq-Gen.
+
+    Parameters
+    ----------
+    seqgen_path : str
+        Path to Seq-Gen executable.
+    seqgen_params : dict
+        Parameter inputs for Seq-Gen.
+    tree : Dendropy.Tree
+    seq_len : int (default: 1000)
+        Lengt of sequences to simulate.
+    gamma_cats : int (default: None)
+        Number of discrete gamma rate categories
+        (continous by default).
+
+    Returns
+    -------
+    matrix, command : tuple (dendropy.CharacterMatrix, str)
+        Simulated dataset and used Seq-Gen command.
+    """
+    s = dendropy.interop.seqgen.SeqGen()
+    s.seqgen_path = seqgen_path
+    s.char_model = 'GTR' if ('general_rates' in seqgen_params) else 'HKY'
+    s.seq_len = seq_len
+    try:
+        s.gamma_shape = seqgen_params['gamma_shape']
+        s.gamma_cats = gamma_cats
+    except KeyError:
+        pass
+    try:
+        s.prop_invar = seqgen_params['prop_invar']
+    except KeyError:
+        pass
+    try:
+        s.state_freqs = seqgen_params['state_freqs']
+    except KeyError:
+        pass
+    if 'general_rates' in seqgen_params:
+        s.general_rates = seqgen_params['general_rates']
+    elif 'ti_tv' in seqgen_params:
+        s.ti_tv = seqgen_params['ti_tv']
+    matrix = s.generate(tree).char_matrices[0]
+    tree_string = tree.as_string('newick', suppress_rooting=True).rstrip()
+    command = ' '.join(s._compose_arguments()) + '\t' + tree_string
+    return (matrix, command)
+
+
+def simulate_multiple_matrices(
+        seqgen_path, p_frame, tree_list, seq_len=1000, gamma_cats=None):
+    """
+    Simulate multiple predictive datasets with Seq-Gen.
+
+    Parameters
+    ----------
+    seqgen_path : str
+        Path to Seq-Gen executable.
+    p_frame : pandas.DataFrame
+        Parameter inputs for Seq-Gen.
+    tree_list : Dendropy.TreeList
+    seq_len : int (default: 1000)
+        Lengt of sequences to simulate.
+    gamma_cats : int (default: None)
+        Number of discrete gamma rate categories
+        (continous by default).
+
+    Returns
+    -------
+    matrix, command : tuple (dendropy.CharacterMatrix, str)
+        Simulated dataset and used Seq-Gen command.
+    """
+    if len(p_frame) == 0:
+        raise ValueError('No parameter values found')
+    if len(p_frame) != len(tree_list):
+        raise ValueError(
+            'Number of parameter values do not match the number of trees.')
+    p_frame['tree'] = tree_list
+    matrices = dendropy.DataSet()
+    seqgen_commands = []
+    for record in _iterrecords(p_frame):
+        tree = record.pop('tree')
+        seqgen_params = get_seqgen_params(record)
+        matrix, command = simulate_matrix(
+            seqgen_path, seqgen_params, tree, seq_len, gamma_cats)
+        matrices.add(matrix)
+        seqgen_commands.append(command)
+    matrices.unify_taxon_namespaces()
+    return (matrices, seqgen_commands)
+
+
+def parse_args(args):
+    parser = argparse.ArgumentParser(
+        description=(
+            'A command-line utility that reads posterior '
+            'output of MrBayes and simulates predictive '
+            'datasets with Seq-Gen.'))
+    parser.add_argument(
+        '-V', '--version', action='version',
+        version='predsim ' + __version__)
+    parser.add_argument(
+        '-l', '--length', type=int, action='store', default=1000,
+        metavar='INT', dest='length',
+        help='sequence lenght (default: 1000)')
+    parser.add_argument(
+        '-g', '--gamma-cats', type=int, action='store', metavar='INT',
+        dest='gamma_cats',
+        help='number of gamma rate categories (default: continuous)')
+    parser.add_argument(
+        '-c', '--commands-file',
+        type=argparse.FileType('w'),
+        dest='commands_file', metavar='PATH',
+        help='path to output file with used Seq-Gen commands')
+    parser.add_argument(
+        '-s', '--skip', type=int, action='store', metavar='INT',
+        dest='skip', default=0, help=(
+            'number of records (trees) to skip at the beginning '
+            'of the sample (default: 0)'))
+    parser.add_argument(
+        '-p', '--seqgen-path',
+        type=str, default='seq-gen',
+        dest='seqgen_path', metavar='PATH',
+        help='path to a Seq-Gen executable (default: "seq-gen")')
+    parser.add_argument(
+        'pfile', type=argparse.FileType('rU'),
+        help='path to a MrBayes p-file')
+    parser.add_argument(
+        'tfile', type=argparse.FileType('rU'),
+        help='path to a MrBayes t-file')
+    parser.add_argument(
+        'outfile', nargs='?', type=argparse.FileType('w'),
+        default=sys.stdout,
+        help='path to output file (default: <stdout>)')
+    return parser.parse_args(args)
+
+
+def main(args=None):
+    if args is None:
+        args = sys.argv[1:]
+    parser = parse_args(args)
+    skiprows = _get_skiprows(parser.skip)
+    p_frame = pandas.read_table(parser.pfile, skiprows=skiprows)
+    tree_list = dendropy.TreeList.get_from_stream(
+        parser.tfile, 'nexus', tree_offset=parser.skip)
+    simulated_matrices, seqgen_commands = simulate_multiple_matrices(
+        parser.seqgen_path,
+        p_frame,
+        tree_list,
+        seq_len=parser.length,
+        gamma_cats=parser.gamma_cats)
+    if parser.commands_file:
+        parser.commands_file.write('\n'.join(seqgen_commands))
+    parser.outfile.write(simulated_matrices.as_string(schema='nexus'))
+
+
+if __name__ == '__main__':  # pragma: no cover
+    main()

requirements.txt

@@ -0,0 +1,7 @@
+pandas>=0.16
+pytest>=2.8
+codecov
+pep8
+pytest-cov
+pytest-pep8
+DendroPy>=4.0

setup.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from setuptools import setup, find_packages
+from os.path import join, dirname
+from io import open
+
+
+setup(
+    name='predsim',
+    version='0.1.0',
+    description=(
+        'Command-line tool for simulating predictive '
+        'datasets from MrBayes\' output.'),
+    long_description=open(
+        join(dirname(__file__), 'README.rst'), encoding='utf-8').read(),
+    packages=find_packages(exclude=['docs', 'tests*']),
+    py_modules=['predsim'],
+    entry_points={
+        'console_scripts': [
+            'predsim = predsim.predsim:main',
+            'predsim.py = predsim.predsim:main']},
+    install_requires=['dendropy>=4.0', 'pandas>=0.16'],
+    author='Markus Englund',
+    author_email='jan.markus.englund@gmail.com',
+    url='https://github.com/jmenglund/predsim',
+    license='MIT',
+    classifiers=[
+        'Development Status :: 5 - Production/Stable',
+        'Intended Audience :: Science/Research',
+        'License :: OSI Approved :: MIT License',
+        'Operating System :: OS Independent',
+        'Programming Language :: Python',
+        'Programming Language :: Python :: 2',
+        'Programming Language :: Python :: 2.7',
+        'Programming Language :: Python :: 3',
+        'Programming Language :: Python :: 3.3',
+        'Programming Language :: Python :: 3.4',
+        'Programming Language :: Python :: 3.5',
+    ],
+    keywords=['simulation', 'Seq-Gen', 'DendroPy'])

test_predsim.py

@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import pytest
+import tempfile
+
+import pandas
+import dendropy
+
+from pandas.util.testing import assert_dict_equal
+
+from predsim import (
+    _get_skiprows,
+    _iterrecords,
+    kappa_to_titv,
+    get_seqgen_params,
+    simulate_matrix,
+    simulate_multiple_matrices,
+    parse_args,
+    main,)
+
+
+SEQGEN_PATH = 'seq-gen'
+
+
+class TestGetSkiprows():
+
+    def test_noskip(self):
+        assert _get_skiprows(0) == [0]
+
+    def test_skip(self):
+        assert _get_skiprows(2) == [0, 2, 3]
+
+
+class TestIterRecords():
+
+    frame = pandas.DataFrame({'a': [1, 2, 3], 'b': [4, 5, 6]})
+
+    def test_iterrecords(self):
+        assert isinstance(list(_iterrecords(self.frame))[0], dict)
+
+
+class TestKappaConversion():
+
+    def test_equal_basefreqs(self):
+        assert kappa_to_titv(1, .25, .25, .25, .25) == 0.5
+
+    def test_unequal_basefreqs(self):
+        assert kappa_to_titv(1, .2, .2, .3, .3) == 0.48
+
+    def test_low_basefreqs_sum(self):
+        assert kappa_to_titv(1, .1, .1, .1, .1) == 0.5
+
+    def test_invalid_basefreqs(self):
+        with pytest.raises(ZeroDivisionError):
+            kappa_to_titv(1, .5, .0, .5, .0)
+
+
+class TestGetSeqGenParamaters():
+
+    d1 = {'pi(A)': 0.25, 'pi(C)': 0.25,	'pi(G)': 0.25,	'pi(T)': 0.25}
+    d2 = {'state_freqs':  '0.25,0.25,0.25,0.25'}
+    d3 = {'pi(A)': 0.25, 'pi(C)': 0.25,	'pi(G)': 0.25}
+
+    def test_equal_basefreqs(self):
+        assert get_seqgen_params(self.d1) == self.d2
+
+    def test_missing_basefreq(self):
+        with pytest.raises(KeyError):
+            get_seqgen_params(self.d3)
+
+
+class TestSingleSimulation():
+
+    tree_string = '((t1:0,t2:0):0,t3:0,t4:0);'
+    tree = dendropy.Tree.get_from_string(tree_string, 'newick')
+
+    rates = {'general_rates': '1,1,1,1,1,1'}
+
+    def test_simulate_empty_params(self):
+        matrix, command = simulate_matrix(SEQGEN_PATH, {}, self.tree)
+        assert len(matrix) == 4
+        assert matrix.sequence_size == 1000
+        assert len(command.split('\t')) == 2
+        assert 'HKY' in command
+
+    def test_gtr(self):
+        matrix, command = simulate_matrix(SEQGEN_PATH, self.rates, self.tree)
+        assert 'GTR' in command
+
+    def test_ti_tv(self):
+        matrix, command = simulate_matrix(
+            SEQGEN_PATH, {'ti_tv': 1}, self.tree)
+        assert 'HKY' in command
+        assert ' -t1 ' in command
+
+    def test_gamma(self):
+        matrix, command = simulate_matrix(
+            SEQGEN_PATH, {'gamma_shape': 2}, self.tree)
+        assert ' -a2 ' in command
+
+
+class TestMultipleSimulations():
+
+    treelist_string = '((t1:0,t2:0):0,t3:0,t4:0);((t1:0,t2:0):0,t3:0,t4:0);'
+    treelist = dendropy.TreeList.get_from_string(treelist_string, 'newick')
+
+    p_frame = pandas.DataFrame({
+        'pi(A)': [0.25, 0.25],
+        'pi(C)': [0.25, 0.25],
+        'pi(G)': [0.25, 0.25],
+        'pi(T)': [0.25, 0.25]})
+
+    def test_empty_input(self):
+        with pytest.raises(ValueError):
+            simulate_multiple_matrices(
+                SEQGEN_PATH, pandas.DataFrame(), dendropy.TreeList())
+
+    def test_two_trees(self):
+        simulate_multiple_matrices(
+            SEQGEN_PATH, self.p_frame, self.treelist)
+
+    def test_parameter_mismatch(self):
+        with pytest.raises(ValueError):
+            simulate_multiple_matrices(
+                SEQGEN_PATH, self.p_frame[:1], self.treelist[:2])
+
+
+class TestArgumentParser():
+
+    def test_parser_help(self):
+        with pytest.raises(SystemExit):
+            parse_args(['-h'])
+
+    def test_parser(self):
+        with tempfile.NamedTemporaryFile() as p_file:
+            with tempfile.NamedTemporaryFile() as t_file:
+                with tempfile.NamedTemporaryFile() as command_file:
+                    parse_args([
+                        '-l100', '-s1', '-g4', '-c', command_file.name,
+                        p_file.name, t_file.name])
+
+
+class TestMain():
+
+    def test_args_help(self):
+        with pytest.raises(SystemExit):
+            main(['-h'])
+
+    def test_noargs(self):
+        with pytest.raises(SystemExit):
+            main()