How to calculate alpha-diversity indice (S.chao1) for each sample and combines with the metadata (Months)?

OSDIAZ · October 7, 2019, 7:49pm

Hi everybody how are you?.
I would be very grateful if you could help me solve this.

I need to calculate Chao1 indixe for each sample and combines with the metadata, specifically with the months.

I already calculate other diversity indeces, and everything was fine, until I try to calculate Chao1.

My script for an example index: SIMPSON:

library(vegan)
library(readxl)
#rarefy con la media de reads por muestra Evaluation of Subsampling-Based Normalization Strategies for Tagged High-Throughput Sequencing Data Sets from Gut Microbiomes - PMC
genero<- read_excel("~/RSTUDIO/Datos_cianobacterias/Cianobacterias_genero_sin3h5c.xlsx")
#columna con las muestrasias
data <- genero
replicates <- as.data.frame(colnames(data)[-1])
colnames(replicates) <- "replicates"
attach(genero)
rwnames <- index
data <- as.matrix(data[,-1])
rownames(data) <- rwnames
data[is.na(data)] <- 0
rarefy(data, 94)
totales <- rowSums(data)
totales <- unname(totales)
mediana <- rowSums(data)
mediana <- unname(totales)
submuestras10 <- rrarefy(data, mediana)
write.csv(submuestras10, file = "~/RSTUDIO/Datos_cianobacterias/Cianobacterias_genero_rarefired_sin3h5c.csv")

#muestras en una fila
genero <- read.csv("~/RSTUDIO/Datos_cianobacterias/Cianobacterias_genero_rarefired_sin3h5c.csv", row.names=1)
data <- submuestras10
#data <- data[,colMeans(data) >=.1]
tot <- rowSums(data) #30 están no rarefied
simpson <- diversity(data, index = "simpson", base = exp(1))
write.csv(simpson, file = "~/RSTUDIO/Datos_cianobacterias/Cianobacterias_genero_rarefired_simpson_sin3h5c.csv")

#boxplot
#agregar el Tx= Control MS o T2DM
bpdata <- read.csv("~/RSTUDIO/Datos_cianobacterias/Metadata-final-sin3h5c.csv", row.names=1)
attach(bpdata)
colvec<- c("yellowgreen","darkturquoise", "coral")
boxplot(simpson~Month, bpdata, col = colvec, ylab = "Simpson", main = "Species Diversity")
library(ggplot2)
library(ggsignif)
res.aov =aov(simpson ~ Month, data = bpdata)
summary(res.aov)#p-value<0.05
TukeyHSD(res.aov)

bpdata$Month <- factor(bpdata$Month,
labels = c("Jul", "Aug", "Sep"))

tumeans <- aggregate(simpson ~ Month, bpdata, mean)
library(ggplot2)
library(ggsignif)
library(tidyverse)
bxplot <- ggplot(bpdata, aes(x=Month, y = simpson, fill =Month )) +
geom_boxplot() +
geom_signif(comparisons = list(c("Jul", "Aug")),
map_signif_level=TRUE,
margin_top = 0.05)+
geom_signif(comparisons = list(c("Aug", "Sep")),
map_signif_level=TRUE,
margin_top = 0.15)+
geom_signif(comparisons = list(c("Jul", "Sep")),
map_signif_level=TRUE,
margin_top = 0.25)+
geom_boxplot(show.legend = FALSE) +
scale_y_continuous(expand = c(0.05,0.05))+
ggtitle("Simpson")+
theme_bw() +
theme(panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(hjust = 0.5, size = 14, family = "Tahoma", face = "bold"),
text=element_text(family = "Tahoma"),
axis.title = element_text(face="bold"),
axis.text.x = element_text(colour="black", size = 11),
axis.text.y = element_text(colour="black", size = 9),
axis.line = element_line(size=0.5, colour = "black"))
scale_fill_brewer(palette = "Accent")
bxplot

My reproducible data:

metadata <- data.frame(data.frame(stringsAsFactors=FALSE,
                                SampleID = c("1A", "1B", "1C", "1D", "1E", "1F", "1G", "1H", "2A",
                                             "2B", "2C", "2D", "2E", "2F",
                                             "2G", "2H", "3A", "3B", "3C", "3D",
                                             "3E", "3F", "3G", "4A", "4B", "4C",
                                             "4D", "4E", "4F", "4G", "4H", "5A",
                                             "5B", "5D", "5E", "5F", "5G", "5H",
                                             "6A", "6B", "6C", "6D", "6E", "6F",
                                             "6G", "6H", "7A", "7B", "7C", "7D",
                                             "7E", "7F", "7G", "7H", "8A", "8B",
                                             "8C", "8D", "8E", "8F", "8G", "8H",
                                             "9A", "9B", "9C", "9D", "9E", "9F",
                                             "9G", "9H", "10A", "10B", "10C", "10D",
                                             "10E", "10F", "10G", "10H", "11A",
                                             "11B", "11C", "11D", "11E", "11F",
                                             "11G", "11H", "12A", "12B", "12C",
                                             "12D", "12E", "12F", "12G", "12H"),
                           SamplingPoint = c("CAJ-1", "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                             "CAJ-2", "CAJ-4", "CAJ-1",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5", "CAJ-2", "CAJ-4", "CAJ-1",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5", "CAJ-2", "CAJ-1", "CAJ-3",
                                             "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-1", "CAJ-3", "CAJ-5", "CAJ-2",
                                             "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                             "CAJ-1", "CAJ-3", "CAJ-5", "CAJ-2",
                                             "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5", "CAJ-1", "CAJ-3", "CAJ-5",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5", "CAJ-2", "CAJ-3", "CAJ-5",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5", "CAJ-2", "CAJ-3", "CAJ-5",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5", "CAJ-2", "CAJ-4", "CAJ-5",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5", "CAJ-2", "CAJ-4", "CAJ-1",
                                             "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                             "CAJ-5"),
                                   Depth = c("80 cm", "80 cm", "Interstitial", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "80 cm",
                                             "80 cm", "Interstitial", "Interstitial",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "Interstitial",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "Interstitial", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "Interstitial", "Interstitial", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "Interstitial",
                                             "Interstitial", "Interstitial", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "Interstitial",
                                             "Interstitial", "Interstitial",
                                             "Interstitial", "Interstitial",
                                             "Interstitial", "80 cm", "80 cm", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "Interstitial", "Interstitial",
                                             "80 cm", "Interstitial", "80 cm",
                                             "80 cm", "80 cm", "80 cm", "80 cm",
                                             "Interstitial", "80 cm", "Interstitial",
                                             "Interstitial", "Interstitial",
                                             "Interstitial", "Interstitial",
                                             "Interstitial"),
                                   Month = c("July", "July", "July", "August", "August", "August",
                                             "September", "September", "July",
                                             "July", "July", "August", "August",
                                             "August", "September", "September",
                                             "July", "July", "July", "August",
                                             "August", "August", "September",
                                             "July", "July", "July", "August",
                                             "August", "August", "September",
                                             "September", "July", "July", "August",
                                             "August", "August", "September",
                                             "September", "July", "July", "July", "August",
                                             "August", "August", "September",
                                             "September", "July", "July", "July",
                                             "August", "August", "September",
                                             "September", "September", "July",
                                             "July", "July", "August", "August",
                                             "September", "September", "September",
                                             "July", "July", "July", "August",
                                             "August", "September", "September",
                                             "September", "July", "July", "July",
                                             "August", "August", "September",
                                             "September", "September", "July", "July",
                                             "July", "August", "August",
                                             "September", "September", "September",
                                             "July", "July", "August", "August",
                                             "August", "September", "September",
                                             "September")
                        )
)
results <- data.frame(data.frame(stringsAsFactors=FALSE,
                                           index = c("1A", "1B", "1C", "1D", "1E", "1F", "1G", "1H",
                                                     "2A", "2B", "2C", "2D",
                                                     "2E", "2F", "2G", "2H", "3A",
                                                     "3B", "3C", "3D", "3E",
                                                     "3F", "3G", "4A", "4B", "4C",
                                                     "4D", "4E", "4F", "4G", "4H",
                                                     "5A", "5B", "5D", "5E",
                                                     "5F", "5G", "5H", "6A", "6B",
                                                     "6C", "6D", "6E", "6F", "6G",
                                                     "6H", "7A", "7B", "7C",
                                                     "7D", "7E", "7F", "7G", "7H",
                                                     "8A", "8B", "8C", "8D", "8E",
                                                     "8F", "8G", "8H", "9A",
                                                     "9B", "9C", "9D", "9E", "9F",
                                                     "9G", "9H", "10A", "10B",
                                                     "10C", "10D", "10E", "10F",
                                                     "10G", "10H", "11A", "11B",
                                                     "11C", "11D", "11E", "11F",
                                                     "11G", "11H", "12A", "12B",
                                                     "12C", "12D", "12E", "12F",
                                                     "12G", "12H"),
                           Un.Caenarcaniphilales = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 2, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 2,
                                                     0, 0, 0, 0),
                          Un.Gastranaerophilales = c(0, 2, 0, 0, 0, 27, 0, 0, 22, 0, 0, 41, 0, 0, 9,
                                                     11, 0, 0, 0, 0, 0, 0, 0,
                                                     203, 0, 0, 28, 0, 0, 2, 15,
                                                     0, 0, 0, 0, 8, 0, 0, 9, 0, 0,
                                                     0, 0, 31, 0, 0, 0, 0, 0, 0,
                                                     2, 0, 0, 0, 19, 8, 0, 24,
                                                     0, 12, 0, 0, 0, 0, 0, 8, 0,
                                                     0, 0, 0, 0, 0, 38, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 14, 0, 0, 0,
                                                     0, 0, 0, 21, 0, 0, 40, 0, 3),
                            Un.Obscuribacterales = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 2, 0,
                                                     0, 0, 0, 0),
                           Un.Vampirovibrionales = c(0, 0, 0, 300, 0, 0, 0, 0, 221, 0, 0, 0, 282, 0,
                                                     0, 7, 232, 0, 0, 95, 94,
                                                     0, 0, 626, 0, 0, 0, 0, 0,
                                                     171, 0, 1139, 0, 0, 0, 0, 57,
                                                     0, 311, 1003, 0, 0, 119, 7,
                                                     99, 0, 422, 0, 0, 0, 0, 0,
                                                     289, 0, 322, 114, 0, 2, 5, 0,
                                                     340, 4, 0, 0, 98, 0, 0, 0,
                                                     0, 9, 121, 152, 0, 0, 0, 105,
                                                     0, 0, 0, 0, 0, 0, 2, 20, 0,
                                                     0, 48, 0, 0, 5, 0, 0, 0,
                                                     15),
                              Un.Melainabacteria = c(0, 0, 0, 0, 273, 470, 0, 0, 363, 110, 169, 0,
                                                     129, 0, 9, 209, 0, 0, 0,
                                                     0, 42, 64, 94, 1050, 154, 0,
                                                     209, 377, 0, 3, 0, 88, 0,
                                                     143, 0, 85, 90, 21, 491, 329,
                                                     45, 0, 215, 384, 0, 154, 21,
                                                     0, 0, 42, 2, 17, 7, 0, 291,
                                                     299, 27, 352, 74, 354, 241,
                                                     2, 0, 0, 0, 749, 0, 116, 0,
                                                     0, 163, 73, 0, 0, 386, 132,
                                                     38, 0, 0, 0, 0, 488, 176,
                                                     0, 0, 0, 0, 71, 40, 273, 0,
                                                     534, 485, 0),
                                 Limnotrichaceae = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 2, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0),
                               Cyanobacteriaceae = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0),
                                  Microcystaceae = c(0, 0, 0, 0, 0, 2, 0, 112, 0, 0, 0, 0, 0, 3, 0, 0,
                                                     0, 0, 0, 0, 0, 2, 0, 7, 0,
                                                     0, 0, 3, 0, 0, 3, 0, 0, 0,
                                                     0, 5, 2, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 55, 0, 0, 0, 0, 0, 2,
                                                     0, 0, 0, 0, 0, 3, 0, 0, 93,
                                                     35, 0, 0, 0, 0, 0, 2, 13,
                                                     396, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 2, 0, 0, 0, 7, 0,
                                                     0, 20, 0, 0, 70, 222),
                                     Nostocaceae = c(0, 0, 19, 0, 1662, 2090, 372, 433, 1000, 769,
                                                     639, 510, 319, 501, 930,
                                                     1876, 316, 613, 0, 0, 1039,
                                                     388, 982, 3272, 895, 243,
                                                     1621, 1094, 43, 291, 2830, 2,
                                                     0, 1006, 2, 1051, 1042, 244,
                                                     1258, 1190, 266, 63, 1016,
                                                     2226, 206, 2014, 183, 0, 27,
                                                     139, 1266, 232, 1036, 102,
                                                     1096, 2439, 39, 2214, 1709,
                                                     2633, 2602, 1403, 2, 5, 46,
                                                     3157, 342, 1031, 407, 5477,
                                                     654, 343, 358, 38, 2265, 824,
                                                     338, 151, 11, 27, 0, 1108,
                                                     1105, 0, 0, 88, 441, 404,
                                                     623, 3218, 172, 2631, 3524,
                                                     1832)
                       )
)

^{Created on 2019-10-07 by the reprex package (v0.3.0)}

This is an example that i found on internet, and i like it, but I could'nt do it

github.com

KasperSkytte/ampvis2/blob/main/R/amp_alphadiv.R

#' Alpha-diversity analysis
#'
#' Calculate alpha-diversity indices for each sample and combines with the metadata.
#'
#' @param data (\emph{required}) Data list as loaded with \code{\link{amp_load}}.
#' @param measure Alpha-diversity measure(s) to be included if not all. A vector of one or more of:
#' \itemize{
#'   \item \code{"observed"}
#'   \item \code{"shannon"}
#'   \item \code{"simpson"}
#'   \item \code{"invsimpson"}
#' }
#' @param richness (\emph{logical}) Also calculate sample richness estimates (Chao1 and ACE) as calculated by \code{\link[vegan]{estimateR}}. (\emph{default:} \code{FALSE})
#' @param rarefy Rarefy species richness to this value before calculating alpha diversity and/or richness. Passed directly as the \code{sample} argument to \code{\link[vegan]{rrarefy}}. (\emph{default:} \code{NULL})
#'
#' @export
#' @importFrom dplyr arrange select
#' @importFrom vegan diversity estimateR
#'
#' @details The alpha-diversity indices are calculated per sample using the vegan function \code{\link[vegan]{diversity}}, where the read abundances are first rarefied using \code{\link[vegan]{rrarefy}} by the size of the \code{rarefy} argument. Refer to the vegan documentation for details about the different indices and how they are calculated. If no measure(s) are chosen, all diversity indices will be returned.

This file has been truncated. show original

Thanks in advance

raytong · October 7, 2019, 11:42pm

Hi @OSDIAZ. I would suggest using estimate_richness from phyloseq package.

library(tidyverse)

metadata <- data.frame(data.frame(stringsAsFactors=FALSE,
                                  SampleID = c("1A", "1B", "1C", "1D", "1E", "1F", "1G", "1H", "2A",
                                               "2B", "2C", "2D", "2E", "2F",
                                               "2G", "2H", "3A", "3B", "3C", "3D",
                                               "3E", "3F", "3G", "4A", "4B", "4C",
                                               "4D", "4E", "4F", "4G", "4H", "5A",
                                               "5B", "5D", "5E", "5F", "5G", "5H",
                                               "6A", "6B", "6C", "6D", "6E", "6F",
                                               "6G", "6H", "7A", "7B", "7C", "7D",
                                               "7E", "7F", "7G", "7H", "8A", "8B",
                                               "8C", "8D", "8E", "8F", "8G", "8H",
                                               "9A", "9B", "9C", "9D", "9E", "9F",
                                               "9G", "9H", "10A", "10B", "10C", "10D",
                                               "10E", "10F", "10G", "10H", "11A",
                                               "11B", "11C", "11D", "11E", "11F",
                                               "11G", "11H", "12A", "12B", "12C",
                                               "12D", "12E", "12F", "12G", "12H"),
                                  SamplingPoint = c("CAJ-1", "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                                    "CAJ-2", "CAJ-4", "CAJ-1",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5", "CAJ-2", "CAJ-4", "CAJ-1",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5", "CAJ-2", "CAJ-1", "CAJ-3",
                                                    "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-1", "CAJ-3", "CAJ-5", "CAJ-2",
                                                    "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                                    "CAJ-1", "CAJ-3", "CAJ-5", "CAJ-2",
                                                    "CAJ-4", "CAJ-1", "CAJ-3", "CAJ-5",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5", "CAJ-1", "CAJ-3", "CAJ-5",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5", "CAJ-2", "CAJ-3", "CAJ-5",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5", "CAJ-2", "CAJ-3", "CAJ-5",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5", "CAJ-2", "CAJ-4", "CAJ-5",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5", "CAJ-2", "CAJ-4", "CAJ-1",
                                                    "CAJ-2", "CAJ-4", "CAJ-1", "CAJ-3",
                                                    "CAJ-5"),
                                  Depth = c("80 cm", "80 cm", "Interstitial", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "80 cm",
                                            "80 cm", "Interstitial", "Interstitial",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "Interstitial",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "Interstitial", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "Interstitial", "Interstitial", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "Interstitial",
                                            "Interstitial", "Interstitial", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "Interstitial",
                                            "Interstitial", "Interstitial",
                                            "Interstitial", "Interstitial",
                                            "Interstitial", "80 cm", "80 cm", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "Interstitial", "Interstitial",
                                            "80 cm", "Interstitial", "80 cm",
                                            "80 cm", "80 cm", "80 cm", "80 cm",
                                            "Interstitial", "80 cm", "Interstitial",
                                            "Interstitial", "Interstitial",
                                            "Interstitial", "Interstitial",
                                            "Interstitial"),
                                  Month = c("July", "July", "July", "August", "August", "August",
                                            "September", "September", "July",
                                            "July", "July", "August", "August",
                                            "August", "September", "September",
                                            "July", "July", "July", "August",
                                            "August", "August", "September",
                                            "July", "July", "July", "August",
                                            "August", "August", "September",
                                            "September", "July", "July", "August",
                                            "August", "August", "September",
                                            "September", "July", "July", "July", "August",
                                            "August", "August", "September",
                                            "September", "July", "July", "July",
                                            "August", "August", "September",
                                            "September", "September", "July",
                                            "July", "July", "August", "August",
                                            "September", "September", "September",
                                            "July", "July", "July", "August",
                                            "August", "September", "September",
                                            "September", "July", "July", "July",
                                            "August", "August", "September",
                                            "September", "September", "July", "July",
                                            "July", "August", "August",
                                            "September", "September", "September",
                                            "July", "July", "August", "August",
                                            "August", "September", "September",
                                            "September")
)
)
results <- data.frame(data.frame(stringsAsFactors=FALSE,
                                 index = c("1A", "1B", "1C", "1D", "1E", "1F", "1G", "1H",
                                           "2A", "2B", "2C", "2D",
                                           "2E", "2F", "2G", "2H", "3A",
                                           "3B", "3C", "3D", "3E",
                                           "3F", "3G", "4A", "4B", "4C",
                                           "4D", "4E", "4F", "4G", "4H",
                                           "5A", "5B", "5D", "5E",
                                           "5F", "5G", "5H", "6A", "6B",
                                           "6C", "6D", "6E", "6F", "6G",
                                           "6H", "7A", "7B", "7C",
                                           "7D", "7E", "7F", "7G", "7H",
                                           "8A", "8B", "8C", "8D", "8E",
                                           "8F", "8G", "8H", "9A",
                                           "9B", "9C", "9D", "9E", "9F",
                                           "9G", "9H", "10A", "10B",
                                           "10C", "10D", "10E", "10F",
                                           "10G", "10H", "11A", "11B",
                                           "11C", "11D", "11E", "11F",
                                           "11G", "11H", "12A", "12B",
                                           "12C", "12D", "12E", "12F",
                                           "12G", "12H"),
                                 Un.Caenarcaniphilales = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 2, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                           0, 0, 0, 0, 0, 0, 0, 0, 2,
                                                           0, 0, 0, 0),
                                 Un.Gastranaerophilales = c(0, 2, 0, 0, 0, 27, 0, 0, 22, 0, 0, 41, 0, 0, 9,
                                                            11, 0, 0, 0, 0, 0, 0, 0,
                                                            203, 0, 0, 28, 0, 0, 2, 15,
                                                            0, 0, 0, 0, 8, 0, 0, 9, 0, 0,
                                                            0, 0, 31, 0, 0, 0, 0, 0, 0,
                                                            2, 0, 0, 0, 19, 8, 0, 24,
                                                            0, 12, 0, 0, 0, 0, 0, 8, 0,
                                                            0, 0, 0, 0, 0, 38, 0, 0, 0,
                                                            0, 0, 0, 0, 0, 14, 0, 0, 0,
                                                            0, 0, 0, 21, 0, 0, 40, 0, 3),
                                 Un.Obscuribacterales = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                          0, 0, 0, 0, 0, 0, 0, 2, 0,
                                                          0, 0, 0, 0),
                                 Un.Vampirovibrionales = c(0, 0, 0, 300, 0, 0, 0, 0, 221, 0, 0, 0, 282, 0,
                                                           0, 7, 232, 0, 0, 95, 94,
                                                           0, 0, 626, 0, 0, 0, 0, 0,
                                                           171, 0, 1139, 0, 0, 0, 0, 57,
                                                           0, 311, 1003, 0, 0, 119, 7,
                                                           99, 0, 422, 0, 0, 0, 0, 0,
                                                           289, 0, 322, 114, 0, 2, 5, 0,
                                                           340, 4, 0, 0, 98, 0, 0, 0,
                                                           0, 9, 121, 152, 0, 0, 0, 105,
                                                           0, 0, 0, 0, 0, 0, 2, 20, 0,
                                                           0, 48, 0, 0, 5, 0, 0, 0,
                                                           15),
                                 Un.Melainabacteria = c(0, 0, 0, 0, 273, 470, 0, 0, 363, 110, 169, 0,
                                                        129, 0, 9, 209, 0, 0, 0,
                                                        0, 42, 64, 94, 1050, 154, 0,
                                                        209, 377, 0, 3, 0, 88, 0,
                                                        143, 0, 85, 90, 21, 491, 329,
                                                        45, 0, 215, 384, 0, 154, 21,
                                                        0, 0, 42, 2, 17, 7, 0, 291,
                                                        299, 27, 352, 74, 354, 241,
                                                        2, 0, 0, 0, 749, 0, 116, 0,
                                                        0, 163, 73, 0, 0, 386, 132,
                                                        38, 0, 0, 0, 0, 488, 176,
                                                        0, 0, 0, 0, 71, 40, 273, 0,
                                                        534, 485, 0),
                                 Limnotrichaceae = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 2, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                     0, 0, 0, 0),
                                 Cyanobacteriaceae = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                       0, 0, 0, 0),
                                 Microcystaceae = c(0, 0, 0, 0, 0, 2, 0, 112, 0, 0, 0, 0, 0, 3, 0, 0,
                                                    0, 0, 0, 0, 0, 2, 0, 7, 0,
                                                    0, 0, 3, 0, 0, 3, 0, 0, 0,
                                                    0, 5, 2, 0, 0, 0, 0, 0, 0,
                                                    0, 0, 55, 0, 0, 0, 0, 0, 2,
                                                    0, 0, 0, 0, 0, 3, 0, 0, 93,
                                                    35, 0, 0, 0, 0, 0, 2, 13,
                                                    396, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                    0, 0, 0, 2, 0, 0, 0, 7, 0,
                                                    0, 20, 0, 0, 70, 222),
                                 Nostocaceae = c(0, 0, 19, 0, 1662, 2090, 372, 433, 1000, 769,
                                                 639, 510, 319, 501, 930,
                                                 1876, 316, 613, 0, 0, 1039,
                                                 388, 982, 3272, 895, 243,
                                                 1621, 1094, 43, 291, 2830, 2,
                                                 0, 1006, 2, 1051, 1042, 244,
                                                 1258, 1190, 266, 63, 1016,
                                                 2226, 206, 2014, 183, 0, 27,
                                                 139, 1266, 232, 1036, 102,
                                                 1096, 2439, 39, 2214, 1709,
                                                 2633, 2602, 1403, 2, 5, 46,
                                                 3157, 342, 1031, 407, 5477,
                                                 654, 343, 358, 38, 2265, 824,
                                                 338, 151, 11, 27, 0, 1108,
                                                 1105, 0, 0, 88, 441, 404,
                                                 623, 3218, 172, 2631, 3524,
                                                 1832)
)
)

results %>%
  as.data.frame() %>%
  column_to_rownames("index") %>%
  as.matrix() %>%
  phyloseq::otu_table(taxa_are_rows = FALSE) %>%
  phyloseq::estimate_richness(measures = "Chao1")
#> Warning in phyloseq::estimate_richness(., measures = "Chao1"): The data you have provided does not have
#> any singletons. This is highly suspicious. Results of richness
#> estimates (for example) are probably unreliable, or wrong, if you have already
#> trimmed low-abundance taxa from the data.
#> 
#> We recommended that you find the un-trimmed data and retry.
#>      Chao1 se.chao1
#> X1A      0      NaN
#> X1B      1        0
#> X1C      1        0
#> X1D      1        0
#> X1E      2        0
#> X1F      4        0
#> X1G      1        0
#> X1H      2        0
#> X2A      4        0
#> X2B      2        0
#> X2C      2        0
#> X2D      3        0
#> X2E      3        0
#> X2F      2        0
#> X2G      3        0
#> X2H      4        0
#> X3A      2        0
#> X3B      1        0
#> X3C      0      NaN
#> X3D      1        0
#> X3E      3        0
#> X3F      3        0
#> X3G      2        0
#> X4A      5        0
#> X4B      2        0
#> X4C      1        0
#> X4D      3        0
#> X4E      3        0
#> X4F      1        0
#> X4G      4        0
#> X4H      3        0
#> X5A      3        0
#> X5B      0      NaN
#> X5D      2        0
#> X5E      1        0
#> X5F      4        0
#> X5G      4        0
#> X5H      2        0
#> X6A      4        0
#> X6B      3        0
#> X6C      2        0
#> X6D      1        0
#> X6E      3        0
#> X6F      4        0
#> X6G      2        0
#> X6H      3        0
#> X7A      3        0
#> X7B      0      NaN
#> X7C      1        0
#> X7D      2        0
#> X7E      3        0
#> X7F      3        0
#> X7G      3        0
#> X7H      1        0
#> X8A      4        0
#> X8B      4        0
#> X8C      2        0
#> X8D      5        0
#> X8E      3        0
#> X8F      3        0
#> X8G      5        0
#> X8H      4        0
#> X9A      1        0
#> X9B      1        0
#> X9C      2        0
#> X9D      3        0
#> X9E      1        0
#> X9F      3        0
#> X9G      2        0
#> X9H      4        0
#> X10A     3        0
#> X10B     3        0
#> X10C     2        0
#> X10D     1        0
#> X10E     2        0
#> X10F     3        0
#> X10G     2        0
#> X10H     1        0
#> X11A     1        0
#> X11B     1        0
#> X11C     0      NaN
#> X11D     3        0
#> X11E     4        0
#> X11F     1        0
#> X11G     0      NaN
#> X11H     1        0
#> X12A     3        0
#> X12B     2        0
#> X12C     4        0
#> X12D     5        0
#> X12E     1        0
#> X12F     3        0
#> X12G     3        0
#> X12H     4        0

^{Created on 2019-10-07 by the reprex package (v0.3.0)}

However, I can't get the meaning of "combines with the metadata".

OSDIAZ · October 8, 2019, 4:24am

Thanks for your help @raytong

I found another method but I'm stuck now.

#this is the sum of my results of every row. (94)
library(vegan)
familia <- matrix(nrow=1,c(8685, 5832, 2736, 7571, 13323, 26618, 4945, 16776, 26950, 18292, 20594, 8516, 10979, 11302, 16243, 25138, 11467, 1934, 3306, 4848, 6322, 6579, 8536, 40754, 22397, 6915, 12979, 12663, 1406, 17175, 45241, 17120, 5075, 11672, 11031, 18511, 7720, 18613, 30095, 30008, 7804, 9210, 9575, 22326, 16090, 53929, 12897, 2878, 4406, 28670, 24283, 16127, 5708, 5912, 24364, 39973, 11233, 26509, 12474, 31516, 18542, 21242, 5768, 11050, 4034, 30936, 24928, 56063, 18179, 15531, 18602, 17354, 4848, 13901, 14953, 15301, 4150, 3866, 4015, 1697, 6952, 16781, 14023, 4846, 8920, 15852, 32939, 13868, 21618, 19606, 7513, 23231, 44961, 44623
))
num_species=specnumber(TR)
chao1=estimateR(TR)[2,]
shannon=diversity(TR,"shannon")
rarecurve(TR)
estimateR(TR)
write.csv(TR, file = "~/RSTUDIO/Datos_cianobacterias/Cianobacterias_familia_TR_sin3h5c.csv")

I don't know how to transpose the final row into column, I need the same lenght as my database in this case months.. (94 DATA )

#boxplot
#this is the database that I already shared with u.
metadata <- read.csv("~/RSTUDIO/Datos_cianobacterias/Metadata-final-sin3h5c.csv", row.names=1)
attach(metadata)
colvec<- c("yellowgreen","darkturquoise", "coral")
boxplot(TR~Month, metadata, col = colvec, ylab = "Chao", main = "Species Diversity")
library(ggplot2)
library(ggsignif)
res.aov =aov(TR ~ Month, data = metadata)
summary(res.aov)#p-value<0.05
TukeyHSD(res.aov)

I don't know if you have any suggest about this.

Thanks

raytong · October 8, 2019, 6:51am

@OSDIAZ. I don't understand your question. Can you describe more. And what is TR?

OSDIAZ · October 10, 2019, 3:25pm

I already solved this issue thank you so much for all ur help.

This is what I did, if you want to know.

library(vegan)
library(readxl)
#rarefy con la media de reads por muestra Evaluation of Subsampling-Based Normalization Strategies for Tagged High-Throughput Sequencing Data Sets from Gut Microbiomes - PMC
familia<- read_excel("~/RSTUDIO/Datos_cianobacterias/Cianobacterias_familia_sin3h5c.xlsx")
#columna con las muestrasias
data <- familia
replicates <- as.data.frame(colnames(data)[-1])
colnames(replicates) <- "replicates"
attach(familia)
rwnames <- index
data <- as.matrix(data[,-1])
estimateR(data)
ch <- estimateR(data)
ch=estimateR(data)[2,]
write.csv(ch, file = "~/RSTUDIO/Datos_cianobacterias/ch-estimater-familia.csv")

#boxplot
#agregar el Tx= Control MS o T2DM
bpdata <- read.csv("~/RSTUDIO/Datos_cianobacterias/Metadata-final-sin3h5c.csv", row.names=1)
attach(bpdata)
colvec<- c("yellowgreen","darkturquoise", "coral")
boxplot(ch~Month, bpdata, col = colvec, ylab = "Chao1", main = "Species Diversity")
library(ggplot2)
library(ggsignif)
res.aov =aov(ch ~ Month, data = bpdata)
summary(res.aov)#p-value<0.05
TukeyHSD(res.aov)

bpdata$Month <- factor(bpdata$Month,
labels = c("Jul", "Aug", "Sep"))

tumeans <- aggregate(ch ~ Month, bpdata, mean)
library(ggplot2)
library(ggsignif)
library(tidyverse)
bxplot <- ggplot(bpdata, aes(x=Month, y = ch, fill =Month )) +
geom_boxplot() +
geom_signif(comparisons = list(c("Jul", "Aug")),
map_signif_level=TRUE,
margin_top = 0.05)+
geom_signif(comparisons = list(c("Aug", "Sep")),
map_signif_level=TRUE,
margin_top = 0.15)+
geom_signif(comparisons = list(c("Jul", "Sep")),
map_signif_level=TRUE,
margin_top = 0.25)+
geom_boxplot(show.legend = FALSE) +
scale_y_continuous(expand = c(0.05,0.05))+
ggtitle("Chao1")+
theme_bw() +
theme(panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(hjust = 0.5, size = 14, family = "Tahoma", face = "bold"),
text=element_text(family = "Tahoma"),
axis.title = element_text(face="bold"),
axis.text.x = element_text(colour="black", size = 11),
axis.text.y = element_text(colour="black", size = 9),
axis.line = element_line(size=0.5, colour = "black"))
scale_fill_brewer(palette = "Accent")
bxplot

Have a good day

system · October 17, 2019, 3:25pm

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.