This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How can I get the Shannon and observed alpha diversity Index value?

Hi I am generating a plot of control and case sample's diversity indices with the following code:

p <- plot_richness(physeq_rarefied, x="type", measures=c("Shannon", "Observed")) + 
  stat_compare_means(method = "wilcox.test") +
  geom_boxplot(aes(fill = type), notch = TRUE, outlier.size = -1) + 
  scale_fill_manual(values = c("hotpink", "skyblue"))+
  labs(x= "Sample types", y= "Alpha Diversity Measure", fill ="Groups",
       title = "Alpha diversity of control and case samples") 
p
p$layers[1] <- NULL
p

I just get two notch plotes, one for case and another for control samples. Now, I want to get the median value for the case and control groups in a table format and export them. How can I do that? Can anyone please help me?

r
• 3,300 views • link
written by dpc dpc

Hello dpc
I used qiime2 for alpha-beta diversity calculeate.

• 0 views • link
written by genetixs • genetixs

1 answer

Hi,

There isn't a way to get directly:

the median value for the case and control groups in a table format and export them.

You need first to get a table, then perform the median for the groups of samples that you're interested in and then export them.

Please see the example below and try to adapt to your data set.

# import package
library("phyloseq") # microbiome analyses: v.1.30.0
library("dplyr") # data manipulation: v.0.8.5

# import test data - phyloseq class obj
data("GlobalPatterns")

# estimate alpha-diversity: table/data.frame
alpha_div <- estimate_richness(physeq = GlobalPatterns, measures = c("Shannon", "Observed"))

## add metadata columns to alpha-diversity data frame 
metadata <- sample_data(object = GlobalPatterns) %>% 
data.frame(.) # convert to a data.frame

# check if all the rownames between metadata and alpha-diversity match to join them blindly
#prints nothing if they match; otherwise prints error
stopifnot( all( rownames(metadata) == rownames(alpha_div) ) )

# add metadata to alpha-div
alpha_div_metadata <- cbind(alpha_div, metadata)

# determine the median by groups: change the group 'SampleType' to the sample_data column 
#name group that you want to use
alpha_div_metadata_median <- alpha_div_metadata %>% 
group_by(SampleType) %>% 
summarise_at(vars(Observed:Shannon), median) %>% 
as.data.frame(.)

# print results
print(alpha_div_metadata_median)

# SampleType Observed  Shannon
# 1              Feces   2648.0 3.640570
# 2         Freshwater   3429.0 3.372561
# 3 Freshwater (creek)   6386.0 3.552736
# 4               Mock   3130.0 4.006375
# 5              Ocean   3427.0 4.483806
# 6 Sediment (estuary)   4581.0 5.461840
# 7               Skin   3214.0 4.849999
# 8               Soil   6964.0 6.576517
# 9             Tongue   2516.5 3.288761   


# save the results as tsv table
write.table(x = alpha_div_metadata_median, file = "alpha_div_metadata_median.tsv", 
            sep = "\t", row.names = FALSE, quote = FALSE)

I hope this helps.

António

• 0 views • link
written by antonioggsousa antonioggsousa

Thanks, Antonio. It is helpful. Can I also draw your attention to my latest query?

dpc

• 0 views • link
written by dpc dpc

Log in to answer this question.

More posts like this