Skip to content
antagomir edited this page Apr 5, 2015 · 41 revisions

Calculate diversity indices

Alternative ways to estimate diversity, richness, and evenness for the samples (NOTE: oligo.matrix must be in absolute scale, not log10!). Examples with simulated data. See read.profiling to use your own data files:

library(microbiome, quietly = TRUE)
data.directory <- system.file("extdata", package = "microbiome")
oligo.matrix.nolog <- read.profiling(level = "oligo", 
			     data.dir = data.directory, log10 = FALSE)
## Reading /home/lei/R/x86_64-pc-linux-gnu-library/3.1/microbiome/extdata/oligoprofile.tab
# NOTE: use absolute-scale data for diversity calculations (no logarithm!)
my.data <- oligo.matrix.nolog  

# This will return a samples x indices table with 
# richness, evenness and diversity collected in one table
div.table <- estimate.diversity(my.data, diversity.index = "shannon")  
## Error in eval(expr, envir, enclos): could not find function "estimate.diversity"
# Estimate richness, evenness, and diversity separately
di <- diversity(my.data, diversity.index = "shannon")
## Error in diversity(my.data, diversity.index = "shannon"): unused argument (diversity.index = "shannon")
ri <- richness(my.data, det.th = NULL)
## Error in eval(expr, envir, enclos): could not find function "richness"
ev <- evenness(my.data, det.th = NULL)
## Error in eval(expr, envir, enclos): could not find function "evenness"

Diversity boxplot

Produce diversity boxplot for your selected sample groups. NOTE: colors and sample groups are specified before function call. To tune y-axis limits, use 'ylim' argument. For other options, see help(diversity.boxplot).

# Define sample groups 
# Alternatively, read metadata from file. See
# https://github.com/microbiome/microbiome/wiki/reading for details
# sample.groups <- metadata$group
sample.groups <- list()
sample.groups$Group1 <- colnames(my.data)[1:3]
sample.groups$Group2 <- colnames(my.data)[4:6]

# Plot diversity boxplots
res <- diversity.boxplot(my.data, sample.groups, diversity.index = "shannon")
## Error in eval(expr, envir, enclos): could not find function "diversity.boxplot"
# The function also returns the sample groups and diversity values 
# used for the plot
sample.groups <- res$sample.groups
div.table <- res$diversity.table 

Writing diversity table into file:

output.dir <- "./"
write.table(div.table,
	file = paste(output.dir, "/DiversityTable.tab", sep=""), sep = "\t")

Phylotype-specific diversity tables

Retrieve phylotypes x samples table which presents diversity within each higher-level taxonomic category for each sample.

library(microbiome)

# Read oligo-level data
data.directory <- system.file("extdata", package = "microbiome")
oligo.matrix.nolog <- read.profiling(level = "oligo", 
			     data.dir = data.directory, log10 = FALSE)  			    	
## Reading /home/lei/R/x86_64-pc-linux-gnu-library/3.1/microbiome/extdata/oligoprofile.tab
# Read mappings between taxonomic levels
phylogeny.info <- GetPhylogeny("HITChip", "filtered")
## Reading /home/lei/R/x86_64-pc-linux-gnu-library/3.1/microbiome/extdata/phylogeny.filtered.tab
# Produce phylotypes vs. samples diversity table
divtab <- diversity.table(oligo.matrix.nolog, phylogeny.info, level.from = "L1", level.to = "oligo", diversity.index = "shannon") 
## Error in eval(expr, envir, enclos): could not find function "diversity.table"

Estimating relative abundancies

# NOTE: estimate relative abundancies for phylotypes based on 
# absolute-scale data for diversity calculations (no logarithm!)
rel <- relative.abundance(my.data, det.th = NULL)
## Warning in relative.abundance(my.data, det.th = NULL): Applying detection
## threshold at 0.8 quantile: 232.026771597465

Abundance table

Calculate abundancies: discretize matrix to form abundancy table of form j, nj where j is number of counts and nj is number of phylotypes with the corresponding counts. This format is often required by richness estimation:

abu <- make.abundancy.table(log10(oligo.matrix.nolog), 
       			det.th = 0, discretization.resolution = 1)

Clone this wiki locally