In austenapigo/lotus: Host Specificity

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Structural Specificity

• structural.specificity: calculates presence-absence host species richness (the number of hosts a symbionts occupies) or abundance-weighted Shannon’s H diversity index (Shannon and Weaver 1948) to quantify symbiont presence and evenness among hosts.
• null.structural: generates a null occupancy-abundance model by randomizing the community matrix and calculating absolute structural specificity per symbiont within each randomized community.
• relative.structural : calculates the deviance in absolute structural specificity from the null expectation per symbiont per host.

Network Specificity

• network.specificity: calculates the presence-absence Resource Range Index or abundance-weighted Paired Difference Index (Poisot et al. 2012) to quantify the 'strength' of host-symbiont interactions by accounting for all potential hosts a symbiont could occupy in a host community.
• null.network: generates a null occupancy-abundance model by randomizing the community matrix and calculating absolute network specificity per symbiont within each randomized community.
• relative.network: calculates the deviance in absolute network specificity from the null expectation per symbiont per host.

Phylogenetic Specificity

• phylogenetic.specificity: calculates as the presence-absence and abundance-weighted Mean Pairwise Phylogenetic Distance (Webb 2000) to quantify symbiont presence and evenness as a function of host phylogenetic breadth. You must supply a phylogenetic tree and the function will calculate a phylogenetic distance matrix for you.
• relative.phylogenetic: calculates the deviance in observed phylogenetic specificity from the null expectation per symbiont per host. Note: null.phylogenetic is combined within deviance.phylogenetic. The picante package conveniently produces null phylogenetic models and calculates standardized effect sizes.

Beta-Specificity

• beta.specificity: calculates the Sørensen (Diserud and Odegaard 2007) or Morisita-Horn (Chao et al. 2008) Multiple-Assemblage Overlap Measure to quantify symbiont interaction consistency to a given host group (species, genus, family, etc.) across space or time.
• null.beta: generates a null occupancy-abundance model by randomizing the community matrix and calculating beta-specificity per symbiont within each randomized community. deviance.beta: calculates and plots the deviance in observed beta-specificity from the null expectation per symbiont per host.
• relative.beta: calculates the deviance in absolute network specificity from the null expectation per symbiont per host.

Important Notes:

• The community matrix for analysis should be organized with hosts populating the rows, while symbionts populate the columns.
• Plant hosts populate the rows, while fungal symbionts populate the columns. Hosts are labeled by their species name with a period and number identifier (e.g., hostA.1) and helps differentiate host samples that are of the same species. This naming scheme is required because host specificity is quantified at the level of host species and not host samples.
• Host specificity is calculated for a symbiont across the entire host community. For example, a symbiont found in given host would be evaluated for its host specificity relative to all hosts that are present in a given dataset.
• More positive values always indicate a narrower symbiont niche and thus higher host specificity. We've negated (multipled by -1) structural and phylogenetic specificity to make this consistent across all metrics.
• For beta-specificity, labels can vary depending on your questions and taxonomic-level of interest. For example, you might be interested in the beta-specificity of symbionts to a taxonomic family of hosts, rather than host species, and could label hosts as Pinaceae.1, Pinaceae.2, Pinaceae.3, etc., for example.
• If your preferred null model is not represented in bipartite::nullmodel you can use any other function to generate randomized communities (e.g., vegan::permatswap) outside of lotus. Make sure your output is formatted as a list of matrices (match output from bipartite::nullmodel) and you can supply this to the randomized.object argument to the null.structural, null.network or null.beta functions. The iterations argument should match the number of iterations in the provided randomized.object.

Common lotus arguments: - abundance.weighted: Logical. TRUE calculates abundance-weighted metrics per symbiont. FALSE calculates presence-absence metrics per symbiont. - model: Character. Specify whether the null expectation should be approximated as a first-(linear) or second-(quadratic) order function. - trim: Logical. TRUE removes symbionts that occupy one host sample. FALSE keeps all symbionts. Note: We think they should be removed because host specificity from an observation of one host species will always result in the highest host specificity value and cannot be relevatized in a meaningful way to a null expectation. - notify: Logical. TRUE prints the current iteration of the for loop. FALSE supresses notifications.

Installation

Install the latest version of lotus from GitHub use:

library(devtools)
devtools::install_github("austenapigo/lotus")

Load the package

library(lotus)
library(ggplot2)
#library(ggpmisc)
#library(picante)
#library(vegan)
#library(bipartite)

Example Data

These examples use the dataset from Small 1976. The study took place in the Mer Bleue peat bog of Ottawa, Canada in 1973. The paper is a preliminary evaluation of the pollination relationships of the major entomophilous plant species of the Mer Bleue. Description taken from (https://iwdb.nceas.ucsb.edu/html/small_1976.html).

Note: This interaction matrix has been modified such that each host species is composed of two host samples. This change was made to maake this dataset compatible for the measurement of beta-specificity.

lotus::comm.matrix # call data

dim(comm.matrix) # check dimensions

comm.matrix # view data frame

Structural Specificity Example

# Calculate uncorrected host specificity (not relavitized to a null model)
hs.object <- structural.specificity(comm.matrix, abundance.weighted = FALSE, trim = TRUE)
head(hs.object)

# Explore data and evaluate relationships between host specificity and symbiont read abundance
ggplot(data = NULL, aes(x = hs.object$Structural.Specificity)) + 
  geom_density(aes(y=..scaled..)) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 12, color = "black"),
    axis.text.y = element_text(size = 12, color = "black"),
    axis.title.x = element_text(size = 12, margin = margin(t = 5, r = 0, b = 0, l = 0)),
    axis.title.y = element_text(size = 12, margin = margin(t = 0, r = 5, b = 0, l = 0)),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    legend.position = "bottom",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line.x = element_blank(),
    axis.line.y = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 12, face = "bold.italic"),
    text = element_text(),
    aspect.ratio = 0.85
  ) +
  labs(x = "Absolute Structural Specificity\n(Host Species Richness)", y = "Proportion of Observations") 

# Make a new data frame with read abundance
read.abund <- as.data.frame(colSums(comm.matrix)) # Calculate read abundances per symbiont
read.abund.trim <- read.abund[rownames(read.abund) %in% hs.object$Symbiont, ] # trim relative to hs.object
structural.object <- data.frame(hs.object, Read.Abundance = read.abund.trim) # make data frame

# Run a correlation test
cor.test(structural.object$Structural.Specificity, log(structural.object$Read.Abundance)) # seems to be significantly negatively correlated suggesting bias - rarer symbionts are more host-specific 

# Visualize occupancy-abundance relationship
ggplot(data = structural.object, aes(y = Structural.Specificity, x = log(Read.Abundance))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red", formula = y ~ x) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 12, color = "black"),
    axis.text.y = element_text(size = 12, color = "black"),
    axis.title.x = element_text(size = 12, margin = margin(t = 5, r = 0, b = 0, l = 0)),
    axis.title.y = element_text(size = 12, margin = margin(t = 0, r = 5, b = 0, l = 0)),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    legend.position = "bottom",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line.x = element_blank(),
    axis.line.y = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 12, face = "bold.italic"),
    text = element_text(),
    aspect.ratio = 0.85
  ) +
  labs(y = "Absolute Structural Specificity\n(-1 * Host Species Richness)", x = "Log Read Abundance Per Symbiont")

# Randomize community matrix to generate a null model for deviance calculations
# This is not the fastest function ever... Setting notify = TRUE will let you know how far along you are. 
null.structural.object <- null.structural(comm.matrix, iterations = 100, abundance.weighted = FALSE, randomization.method = "shuffle.web", trim = TRUE, notify = TRUE)
head(null.structural.object)

# Calculate and plot the deviance of observed host specificity from the null boundary and get averages per host sample
structural.dev <- relative.structural(comm.matrix, randomized = null.structural.object, abundance.weighted = FALSE, trim = TRUE, notify = TRUE)
head(structural.dev) # View data frame of output

ggplot(data = structural.dev, aes(y = Mean.Standardized.Effect.Size, x = Host.Sample)) +
  geom_pointrange(aes(ymin = Mean.Standardized.Effect.Size - Standard.Error.of.Mean.SES, ymax = Mean.Standardized.Effect.Size + Standard.Error.of.Mean.SES)) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 12, color = "black", angle = 90),
    axis.text.y = element_text(size = 12, color = "black"),
    axis.title.x = element_text(size = 12, margin = margin(t = 5, r = 0, b = 0, l = 0)),
    axis.title.y = element_text(size = 12, margin = margin(t = 0, r = 5, b = 0, l = 0)),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    legend.position = "bottom",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line.x = element_blank(),
    axis.line.y = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 12, face = "bold.italic"),
    text = element_text(),
    aspect.ratio = 0.85
  ) +
  labs(y = "Relative Structural Specificity\n(-1 * Host Species Richness)", x = "Host Sample")

You can read more about each lotus function with the help function

help("structural.specificity")
help("null.structural")
help("relative.structural")

austenapigo/lotus documentation built on Aug. 8, 2021, 10:54 a.m.