# Para Loles :bird:
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.
:::info
`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.
```r=1
#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](https://github.com/laurichi13/PYM) del estudio puedes crear el objeto phyloseq, y también tienes el archivo `.csv` con la [metadata](https://github.com/laurichi13/PYM/blob/main/Data/metadata.csv).
Cualquier duda me dices porque a lo mejor no tengo todo colgado en el github :smile: ...
Para obtener una tabla con todo lo que puedes correlacionar, puedes usar este script (modificando adecuadamente)
```r=46
#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")
```
:::info
:warning: 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
```r=1
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
```r=1
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
```r=1
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](https://microbiome.github.io/tutorials/Heatmap.html) pipeline, a partir de "Cross-correlating data sets"
Puedes fijarte en [este](https://github.com/laurichi13/NRCD_analysis/blob/main/R%20Scripts/complex_heatmap.R) 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](https://microbiome.github.io/OMA/docs/devel/) tutorial es nuestra pequeña Biblia.
- El paquete [mia](https://microbiome.github.io/mia/)
- El paquete [phyloseq](https://joey711.github.io/phyloseq/)
- El paquete [microbiome](https://microbiome.github.io/tutorials/).
- 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](https://docs.qiime2.org/2024.10/tutorials/index.html)
Sign in with Wallet
Connect another wallet