How can I get the Shannon and observed alpha diversity Index value?
1
0
Entering edit mode
10 months ago
dpc ▴ 210

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 • 554 views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
3
Entering edit mode
10 months ago
antonioggsousa ★ 2.4k

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

ADD COMMENT
0
Entering edit mode

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

dpc

ADD REPLY

Login before adding your answer.

Traffic: 1749 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6