Visualization overhaul (2 of 2)

[polyatail]

The PR is a substantial overhaul to how Samples and their associated Classifications are analyzed in a notebook environment.

The following issues were addressed:

• Add PCoA as a plotting option #74 Scatter plots of distances between samples (calculated using a distance metric, e.g., unifrac) can now be visualized with principal component analysis, via two algorithms: deterministic pcoa, or smacof, an iterative optimization. This is exposed to the user via plot_mds.
• Make resource creation or lookup easier #130 Enabled easier lookup of Tags resources only, rather than a change affecting all models.
• Explore static image rendering options for notebooks #145 To facilitate rendering of report PDFs (via HTML and weasyprint), setting an environment variable will enable the onecodex renderer in altair, which produces high-res PNGs.
• Add clustering to sample vs. abundance heatmaps #149 In plot_heatmap, the samples and taxa are ordered according to euclidean distance.
• Allow specification of clustering algorithm in heatmaps, distance plots #150 When clusters are generated, in plot_heatmap or in plot_distance, the user can now select the clustering algorithm, e.g., single- or average-linkage.
• Adopt new docstring conventions and update docstrings #152 New methods and most older methods have updated Numpy-style docstrings now.
• Slices of ResourceList objects should be ResourceList objects #158 Fixed bug where ResourceList slices were returned as lists, when they should have been ResourceLists.
• Return lists of Classifications, Analyses, and Samples as SampleCollection #159 Lists of Samples or Classifications are now returned wrapped in a new SampleCollection object which exposes lots of additional analyses methods that operate on the objects contained within.
• __setattr__ should be aware of ResourceList types #160 Fixed bug where properties of a model which were ResourceList objects were not being saved.

New functionality has been added as part of the SampleCollection object, some examples of which are outlined below.

Retrieving a project and generating plots
All plotting functions are now available as methods of the new SampleCollection object.

>>> project1 = ocx.Projects.get('d53ad03b010542e3')
>>> samples1 = ocx.Samples.where(project=project1.id, public=True, limit=10)
>>> type(samples1)
onecodex.models.SampleCollection
>>> samples1.plot_pca(color='Bacteroides')

Obtaining a DataFrame OTU table
We introduce a new subclassed ResultsDataFrame which is returned when the user requests an OTU table from a SampleCollection. This object transfers taxonomy and metadata along with it, which is automatically subset to only ever contain the taxonomy and metadata of taxa present in the ResultsDataFrame and metadata belonging to Samples present.

results = samples1.results(rank='phylum', normalize=False)

# take the fourth root and plot a heatmap
(results ** 0.25).ocx.plot_heatmap()

# restrict to taxa with 100 reads or more and run PCA
results.loc[:, results.sum() > 100].ocx.plot_pca(color='Bacteroidetes')

# do PCoA with weighted Unifrac
results.ocx.plot_mds(method='pcoa', metric='weighted_unifrac', color='Bacteroidetes')

Metadata and taxonomy can be accessed two ways:

samples1.metadata == samples1.results().ocx.metadata
samples1.taxonomy == samples1.results().ocx.taxonomy

The values in the ResultsDataFrame can be manipulated as usual. Analysis methods which are available to SampleCollection are also available in the ResultsDataFrame.ocx namespace via pd.api.extensions.

Manipulating contents of SampleCollection
The Samples in a SampleCollection can be manipulated as if they were a list. A user might notice that one sample is "bad" and want to exclude it from their analyses. Or they might want to combine multiple projects together, or analyze their own project in the context of a larger public dataset.

project2 = ocx.Projects.get('be33784611fe4c19')
samples2 = ocx.Samples.where(project=project2.id, public=True, limit=10)

# analyze these two projects together
combined = samples1 + samples2

# exclude this sample
del combined[4]

combined.plot_pca(size='Firmicutes')

# get a combined OTU table and write to a CSV file
combined.results(rank='genus', normalize=True).to_csv('otu_table.csv')

# and in long format
combined.results(table_format='long').to_csv('otu_table.long.csv')

[coveralls]

Coverage decreased (-0.6%) to 84.078% when pulling 8c99d6d on roo-viz-part2 into 6473770 on master.

[bovee]

I really like the direction this is going and I think this is a useful reorganization you have here. There are some big changes with some new parts though so more docstrings, even if they're on Mixins or components that aren't otherwise exposed, would be really helpful to understand the overall architecture. Also, I have some comments/concerns about how we're bumping some of our dependencies.

[bovee]

I'm fairly reluctant to pin with a >=; we should probably pin to specific minor-release versions of packages, e.g. using ~=, to prevent unexpected downstream breakages from taking out this library.

[audy]

Here are some comments/suggestions:

Separate normalization steps from stats

I feel pretty strongly that preprocessing/normalization steps should remain separate from statistics steps. This would read nicer and cut down on repetition. The user only has to transform their data once, and then pass the transformed object to multiple stats functions.

# perform normalization

processed_analysis = normalize(agglomerate(analysis_object, rank='Phylum'))

# do some stats
alpha = alpha_diversity(processed_analysis)

beta = beta_diversity(processed_analysis)

API Design

Since the entire API is changing I think it's worthwhile to evaluate a few possible API designs

Proposal A: method chaining. This would be pretty familiar to anyone used to phyloseq (or d3). I also think it reads well and is easy to type:

experiment.agglomerate('Phylum').normalize().beta_diversity()

Proposal B: IDK what this is called but it's familiar to scikit-learn which is an API for processing DataFrames and doing stats/machine learning:

transformer = ReadCountNormalizer(method="proportion")

summary_statistic = ShannonDiversityIndex()

stats = summary_statistic(transformer.transform(results_data_frame))

This API design can be made a little bit more readable using something like Pipeline

# straight up ripoff of scikit-learn
my_analysis = Pipeline([
  ('transform', ReadCountNormalizer()),
  ('top5', SelectTopTaxa(n=5)),
  ('alpha_diversity', ShannonDiversityIndex())
])

my_analysis.transform(results_data_frame)

This design also lets you store artifacts from the transformations / stats functions (e.g., top5.kept_taxa_names)

It would even be possible to have multiple APIs. Again these are just suggestions but I think it's worthwhile. Also, it's possible to implement these APIs on top of what you've built already (and have multiple ones at the same time).

Better naming

"ResultsDataFrame" is kind of ambiguous? I'm bad at naming things though. I am thinking that there are going to be multiple things named "results" in users' notebooks and it might get confusing. The results function is also ambiguous. I think this should be renamed to transform maybe.

Guessing normalization

Even though the odds of a false positive are probably nil, a good alternative would be to store a bool on ResultsDataFrame and have any normalization function turn it to True.

onecodex/onecodex/helpers.py

Lines 107 to 111 in 4e7a769

	def _guess_normalized(self):
	# it's possible that the _results df has already been normalized, which can cause some
	# methods to fail. we must guess whether this is the case
	return bool((self._results.sum(axis=1).round(4) == 1.0).all())

[audy]

(playing around with the lib now. I'll give you some thoughts on my first impressions)

onecodex/onecodex/models/__init__.py

Line 238 in 4e7a769

def primary_classifications(self):

I think @property doesn't work well here because it's really a function that calls an API (and changes the objects state). As a user, I expect that this will return immediately when it's really doing 900 GET requests.

I think it would be better to drop the internal caching mechanisms and just use something external and explicit like https://github.com/reclosedev/requests-cache which automatically caches all HTTP requests and stores them in an sqlite3 database. As a user, I am unaware of when things are being pulled from the API or from the cache. This way, I am in control of the cache. We could make it simpler for the user (like onecodex.enable_cache()).

[polyatail]

IIRC primary_classifications doesn't actually trigger any API calls, but I get what you're saying.

[audy]

SampleCollection.primary_classification checks self._cached and if the result is False, dispatches _classification_fetch which will trigger an API request at least in the case of the model being Classifications.

[audy]

Updated SampleCollections.primary_classifications to only cache the results after all classifications have been fetched to prevent a potential gotcha.

[boydgreenfield]

Minor @polyatail – but this would probably have been a good PR to either rebase or squash before merging into master (commit history starts fairly clean, but also includes a bunch of "More changes..." style commits).