Generate random (gaussian or realistic) distance matrices

.gitignore

@@ -1,3 +1,5 @@
+bayesian-regression/
+
 # OS generated files
 .DS_Store
 .DS_Store?
@@ -53,4 +55,4 @@ pip-log.txt
 # Generated PCoA output of algorithms
 out/*
 
-.env
+.env

environment.yml

@@ -5,13 +5,21 @@ channels:
 - defaults
 - bioconda
 - biocore
+- anaconda
 dependencies:
-- scikit-bio=0.4.2
-- click=6.7
-- python=3
+- scikit-bio
+- click=6.7 # install from anaconda channel
+- python=3.5
 - numpy
 - scikit-learn=0.19.1
 - scipy
 - pandas
 - jupyter
 - matplotlib
+- gneiss
+- pytz
+- tensorflow
+- pip
+# also need to install bayesian-regression:
+- pip:
+  - bayesian_regression-0.1-py3-none-any.whl

scripts/randdm

@@ -0,0 +1,115 @@
+#!/usr/bin/env python
+
+import os
+
+import click
+import numpy as np
+from skbio import DistanceMatrix
+from skbio.stats.distance import randdm
+from bayesian_regression.util.generators import band_table, block_table
+from skbio.diversity import beta_diversity
+
+
+@click.command()
+@click.option('--structure', 'structure', type=str,
+              help='Optional. Generate structure in the random data that '
+                   'models distributions of microbial communities. '
+                   'Options: "band" or "block". If not specified, then '
+                   'draws from gaussian distribution to generate distance '
+                   'matrices.', default=None)
+@click.option('--dim', 'dimensions', type=int, multiple=True,
+              help='Dimension of the input matrix/matrices to generate')
+@click.option('--seed', 'seed', type=int,
+              help='Random number generator seed value.')
+@click.option('--sub', 'subsample_dims', type=str, multiple=True,
+              help='Generate subsampled distance matrix with given'
+                   'dimension (must be smaller than original dimension)'
+                   ' for each randomly generated distance matrix.'
+                   ' Can specify as int or percentage')
+@click.option('--overwrite', 'overwrite', type=bool, default=False,
+              help='Overwrite output directory if it already exists. Do not '
+                   'overwrite, by default.')
+@click.argument('output_dir', type=click.Path(), default=None)
+def generate(dimensions, output_dir, seed, subsample_dims, structure,
+             overwrite):
+    """
+    Generate random distance matrix
+    """
+
+    if seed is not None:
+        np.random.seed(seed)
+
+    if (not overwrite and os.path.exists(output_dir) and not click.confirm(
+            'The output directory %s exists. '
+            'Do you want to overwrite?' % output_dir)):
+        click.echo('Aborting.')
+        return
+
+    if not os.path.exists(output_dir):
+        os.makedirs(output_dir)
+
+    for dim in dimensions:
+        outpath = os.path.join(output_dir, 'randdm-{}.txt'.format(dim))
+
+        click.echo('Generating random distance matrix: %s' % outpath)
+
+        # Generate random distance matrix
+        if structure is not None:
+            click.echo('Generating synthetic OTU table with "{}"-like '
+                       'microbial distribution...'.format(structure))
+            if structure == 'band':
+                biom_table = band_table(dim, dim, seed=seed)[0]
+            elif structure == 'block':
+                biom_table = block_table(dim, dim, seed=seed)[0]
+            else:
+                raise ValueError('Invalid value for --structure parameter.')
+
+            otu_table = biom_table.matrix_data.todense()
+
+            click.echo('Generating distance matrix from OTU table...')
+            distance_matrix = beta_diversity('braycurtis', otu_table,
+                                             validate=False)
+
+        else:
+            click.echo('Generating random distance matrix '
+                       '(gaussian distribution)')
+            distance_matrix = randdm(dim, constructor=DistanceMatrix)
+
+        # Serialize distance matrix
+        distance_matrix.write(outpath)
+
+        # Subsampling
+        for subsample_dim in subsample_dims:
+            # Parse parameter values into integers
+            if '%' in subsample_dim:
+                percent = float(subsample_dim[:-1]) / 100
+                subsample_dim = int(percent * dim)
+            else:
+                subsample_dim = int(subsample_dim)
+
+            if subsample_dim >= dim:
+                raise ValueError('Subsample dimension %d must be smaller than'
+                                 'original matrix dimension %d' % (
+                                     subsample_dim,
+                                     dim))
+            subsampled_outpath = os.path.join(output_dir,
+                                              'randdm-{}-sub-{}.txt'
+                                              .format(dim, subsample_dim))
+
+            click.echo(
+                'Subsampling original randomly generated distance matrix '
+                'with subsample dimension %d: %s' % (subsample_dim,
+                                                     subsampled_outpath))
+
+            ids = distance_matrix.ids
+            # Subsample without replacement
+            subsampled_ids = np.random.choice(ids, subsample_dim,
+                                              replace=False)
+            subsampled_matrix = distance_matrix.filter(subsampled_ids)
+            subsampled_matrix.write(subsampled_outpath)
+
+    click.echo('Done.')
+
+
+if __name__ == '__main__':
+    generate()

setup.py

@@ -44,6 +44,7 @@
       url='http://github.com/biocore/mds-approximations/',
       test_suite='nose.collector',
       packages=find_packages(),
-      scripts=['scripts/mdsa', 'scripts/kruskal', 'scripts/procrustes'],
+      scripts=['scripts/mdsa', 'scripts/kruskal', 'scripts/procrustes',
+               'scripts/randdm'],
       extras_require={'test': ["nose", "pep8", "flake8"]},
       classifiers=classifiers)