Para Loles

Los datos de microbioma los trabajamos con un paquete de R que se llama phyloseq es un paquete que está diseñado para trabajar con datos de micro y poder tener en un mismo "objeto" los datos de abundancia, taxonomía y metadatos.

Phyloseq es un paquete de R que se usa mucho para análisis de microbioma. Un objeto Phyloseq está formado por tres elementos necesarios:
1. La tabla de taxonomía.
2. La tabla de metadata.
3. La tabla de abundancia. Y además…
4. La información filogenética. No es fundamental, pero sin ella no podemos calcular la distancia Unifrac.

Se crea facilmente a partir de la tabla de resultados generada por herramientas de análisis como qiime2, kraken o metaphlan. Para 16S rRNA vienen de qiime2 y puedes usar este script para importar los datos (Con las modificaciones pertinetes). Si ves hago una "limpieza" de la tabla de taxonomia reordenando los taxones para que se lean más facilmente.

#Load libraries
library(tidyverse)
library(qiime2R)
library(phyloseq)
library(microbiome)

#Read data from qiime2####

ASV<-read_qza("Data/table.qza")
#Taxonomy
taxonomyq <-read_qza("Data/taxonomy.qza")

#Transforme table into phyloseq-readable table
taxtableq <- taxonomyq$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))
taxtableq$Kingdom <- gsub("d__", "",taxtableq$Kingdom)
taxtableq$Phylum <- gsub("p__", "",taxtableq$Phylum)
taxtableq$Class <- gsub("c__", "",taxtableq$Class)
taxtableq$Order <- gsub("o__", "",taxtableq$Order)
taxtableq$Family <- gsub("f__", "",taxtableq$Family)
taxtableq$Genus <- gsub("g__", "",taxtableq$Genus)
taxtableq$Species <- gsub("s__", "",taxtableq$Species)

#Load tree
tree<-read_qza("Data/rooted-tree.qza")

#Load Metadata
metadata<-read_csv("Data/metadata_summary.csv")

#Clean-up taxa table
taxtableq <- taxtableq[grep("Unassigned|Eukaryota", taxtableq$Kingdom, invert=T),]
taxtableq <- taxtableq[grep("Chloroplast| Mitochondria|unidentified|uncultured", taxtableq$Family, invert=T),]
taxtableq <- taxtableq[grep("unidentified|uncultured", taxtableq$Genus, invert=T),]

#Create phyloseq object 
OTUs <- otu_table(ASV$data, taxa_are_rows = T) 
tree <- phy_tree(tree$data) 
TAXq <- tax_table(as.data.frame(taxtableq) %>% select_("-Confidence") %>% column_to_rownames("Feature.ID") %>% as.matrix()) #moving the taxonomy to the way phyloseq wants it
sample_metadata <- sample_data(metadata %>% as.data.frame()) 
sample_names(sample_metadata) <- paste(metadata$SampleID)
physeq<- merge_phyloseq(OTUs, tree, TAXq,sample_metadata)


#Prune taxa at genus level and transform data for CoDA  
physeq <- tax_glom(physeq, taxrank = "Genus")
physeq.rel <- transform(physeq, transform="compositional")

Datos de PYM

Con los datos que estan en el github del estudio puedes crear el objeto phyloseq, y también tienes el archivo .csv con la metadata.

Cualquier duda me dices porque a lo mejor no tengo todo colgado en el github

…

Para obtener una tabla con todo lo que puedes correlacionar, puedes usar este script (modificando adecuadamente)

#Write data for sharing
data <- psmelt(physeq) 
ASV_tab <- data %>% select(SampleID,Abundance,Genus) 
genus <- pivot_wider(ASV_tab, id_cols=SampleID, names_from=Genus, values_from=Abundance, values_fn = mean)
write_csv(genus,"genus.csv")

Dependiendo de lo que vayas a hacer tenemos que hablar de filtrados y transformaciones. Alfa diversidad con el objeto sin transformar (physeq) beta diversidad y tablas de abundancias con el physeq.rel.

Para filtrar las muestras con dietas falsas usar este comando:
En vez de Diagnosis sería SampleID = NUMERO

#Remove samples with no diagnosis - Some extra samples were inlcuded by mistake-  
physeq <- subset_samples(physeq, Diagnosis!="NA")

Ejemplos de cáculos de diversidad y abundancias:

Porcentaje de Phylum más abundantes (Top 5 phylum)

Esto lo puedes hacer con género y seleccionando los 20 más abundantes

library(tidyverse)
#Porcentaje Phylum 
Phyl <- tax_glom(pseq.compositional, taxrank = "Phylum") 
tb_phyl <- psmelt(Phyl) 
tbphylum <- tb_phyl %>% group_by(Phylum) %>% summarise(M=median(Abundance)) %>% ungroup() %>% unique() %>% top_n(M,n=5)
tbphylum <- tbphylum %>% inner_join(tb_phyl)
Plotf <- ggplot(tbphylum, aes(SampleID, Abundance ,fill=Phylum)) + geom_col(position="fill") + scale_fill_manual(values=paleta5)
Plotf <- Plotf + theme_minimal() + xlab("SampleID") + ylab("Relative Abundance") +
  theme(axis.text.x = element_text(color = "black", size = 14, angle = 45, hjust = .5, vjust = .5, face = "plain")) +
  theme(axis.text.y = element_text(color = "black", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain")) + 
  theme(legend.title=element_text(color= "black", size=12), legend.text=element_text(size=11, face="italic"))
Plotf <- Plotf +  ggtitle ("") + theme(plot.title = element_text(family="Arial", size=rel(2), vjust=2,face="plain", color="black", lineheight=1.5))  
Plotf

Alfa diversidad

plot_richness(physeq) #paquete phyloseq
rich <- richness(physeq) #paquete microbiome
plot_richness(physeq, color = "Cluster", x = "Cluster", 
              measures = c("Chao1", "Simpson", "Shannon")) + 
              geom_boxplot(aes(fill = Cluster), alpha=.7) + 
              scale_color_manual(values = c("#B2ABD2", "#b2df8a")) + 
              scale_fill_manual(values = c("#B2ABD2", "#b2df8a")) + 
              theme_bw() # phyloseq

Beta diversidad

library(microbiome)
pseq.compositional <- transform(physeq, "compositional")
#Distancia Unifrac
ordinated_taxa_unifrac = ordinate(physeq.composiitonal, method="MDS", distance="unifrac",weighted=TRUE)
plot_ordination(physeq.compositional, ordinated_taxa_unifrac, color="Cluster", title = "Unifrac MDS Analysis") + theme_bw()

Cálculo de correlaciones

Utiliza este pipeline, a partir de "Cross-correlating data sets"

Puedes fijarte en este script mío, pero lo suyo es que sigas las instrucciones que te paso arriba.

Otra info:

• Tutoriales de análisis de mcirobioma útiles: (A partir de la tabla de abundancias)
  • Este tutorial es nuestra pequeña Biblia.
  • El paquete mia
  • El paquete phyloseq
  • El paquete microbiome.
• Para generar la tabla de abundancias desde el raw data 16SrRNA usamos Qiime2, esto si quiere slo vemos más adelante para no crear agobios…
  • tutoriales