r diversity.r

Posted 2021-05-24

get_div <- function(i, physeq){
    # Subsample
    rarefied_physeq <- rarefy_even_depth(physeq, sample.size = min_lib, verbose = FALSE, replace = TRUE)

    # Calculate observed richness for that group and store value in a df column
    richness <- estimate_richness(rarefied_physeq, measures = "Observed")[ ,1]
    # Calculate Simpson's E for that group and store value in a df column
    simpson <- (estimate_richness(rarefied_physeq, measures = "InvSimpson")[ ,1])
    simpson_e <- (estimate_richness(rarefied_physeq, measures = "InvSimpson")[ ,1]/richness)
    shannon <- estimate_richness(rarefied_physeq, measures = "Shannon")[ ,1]
    pielous <- shannon/log(richness)
    df <- data.frame(label = sample_names(rarefied_physeq), iter = rep(i, nsamples(rarefied_physeq)), richness = richness, invsimpson = simpson, shannon = shannon, pielou = pielous, simpson_e = simpson_e)
    return(df)
}

library(phyloseq)
library(pbmcapply)
trials <- 1000
physeq_filt_prev <- prune_taxa(taxa_sums(physeq_filt_prev) > 0, physeq_filt_prev)
min_lib <-   min(sample_sums(physeq_filt_prev))
d16S_alpha <- pbmclapply(1:trials, get_div, physeq = physeq_filt_prev, mc.cores = 24)