The `subset_ord_plot` function is a "convenience function" intended to make it easier to retrieve a plot-derived `data.frame` with a subset of points according to a `threshold` and `method`. The meaning of the `threshold` depends upon the `method`.

Load the necessary packages and data.
```{r message=FALSE, warning=FALSE}
library("phyloseq"); packageVersion("phyloseq")
library("ggplot2"); packageVersion("ggplot2")
data(GlobalPatterns)
```

ggplot2 package theme set. See [the ggplot2 online documentation](http://docs.ggplot2.org/current/) for further help.

```{r ggplot2-theme-set}
theme_set(theme_bw())
```


Some subsetting and light massaging of the `GlobalPatterns` dataset.

Clean zeros.
```{r}
GP <- GlobalPatterns
GP <- prune_species(taxa_sums(GP)>0, GP)
```

Add `"human"` variable to `GP` to indicate human-associated samples.
```{r}
sample_data(GP)$human <- get_variable(GP, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
```

Subset to just Bacteroidetes phylum. 
```{r}
GP <- subset_taxa(GP, Phylum=="Bacteroidetes")
```

Perform a correspondence analysis and then create some plots to demonstrate using `subset_ord_plot`. We want to make species topo with a subset of points layered. Start by performing the correspondence analysis.

```{r}
gpca <- ordinate(GP, "CCA")
```

Now make a basic plot of just the species points in the correspondence analysis.
```{r}
p1 = plot_ordination(GP, gpca, "species", color="Class")
p1
```

Re-draw this as topo without points, and facet
```{r}
p0 = ggplot(p1$data, p1$mapping) + geom_density2d() + facet_wrap(~Class)
p0
```

Re-draw this but include points
```{r}
p1 = p1 + geom_density2d() + facet_wrap(~Class)
p1
```

Add a layer of a subset of species-points that are furthest from origin.
```{r}
p0 + geom_point(data=subset_ord_plot(p1, 0.7, "square"), size=1)
p0 + geom_point(data=subset_ord_plot(p1, 0.7, "farthest"), size=1)
p0 + geom_point(data=subset_ord_plot(p1, 0.7, "radial"), size=1)
```

Here is what the data retreived by `subset_ord_plot` actually looks like
```{r}
head(subset_ord_plot(p1, 0.7, "radial"))
```

