Try   HackMD

Análisis de resultados del gen 16S rRNA con R

Dra. E. Ernestina Godoy Lozano
elizabeth.godoy@insp.mx
tinagodoy@gmail.com

Temario

Crear el espacio de trabajo

Ejercicio 1. Conectarnos al servidor

ssh -X alumnoX@132.248.32.180

Ejercicio 2: Descomprimir Rproject

Ir al directorio EMetagen

cd EMetagen

Descomprimir modulo de trabajo

tar -xzvf /Share/EMetagen/16S_analisis.tar.gz

Ejercicio 3: Conectarnos al R Server

  1. Abrir un navegador web

  2. Ir a la siguiente direccion (link/liga):
    http://132.248.32.180:8787

  3. Ingresar usuario (Username) y contraseña (Password)

    Image Not Showing Possible Reasons
    • The image was uploaded to a note which you don't have access to
    • The note which the image was originally uploaded to has been deleted
    Learn More →

  4. ¡Perfecto! Haz ingresado al R Studio Server!!!

Análisis en R con RStudio

Configuración del espacio de trabajo

Ejercicio 4: Posicionarnos en el directorio de trabajo




getwd()
setwd("~/EMetagen/16S_analisis")
dir()

Ejercicio 5: Cargar librerias de trabajo

Usaremos la función library()










library(gplots)
library(ggplot2)
library(vegan)
library(reshape2)
library(metagenomeSeq)
library(RColorBrewer)
library(phyloseq)
library(MASS)
library(gtools)

Cargar funciones extras

Usaremos la función source()




source("scripts/pairwise.adonis.R")
source("scripts/abund_relativa.R")
source("scripts/taxonomic_levels.R")

Matrices de abundancia

Ejercicio 6: Modificación de matriz de abundancia

Niveles taxonómicos

Diapositivas: 15-17

  1. Cargar datos a espacio de trabajo






matrix_abund <- read.delim("data/TablaFinal.txt", header=TRUE)
View(matrix_abund)
dim(matrix_abund)
tabla <- as.matrix(data.frame(matrix_abund[,2:16]))
rownames(tabla) <- matrix_abund$taxonomy
View(tabla)
  1. Escribir matriz de abundancias absolutas modificada
exportMat(tabla, file ="data/integrated_matrix.txt")
  1. Añadir los nombres de los taxa a un objeto


asignacion_taxonomica <- data.frame(completo=rownames(tabla))
View(asignacion_taxonomica)
  1. Asignar los taxa a su nivel taxonomico correspondiente
    Usaremos la función sapply() para aplicar una función a todo un dataFrame







asignacion_taxonomica$Domain <- sapply(as.character(asignacion_taxonomica$completo), function(x){strsplit(x, ";")[[1]][1]})
asignacion_taxonomica$Phylum <- sapply(as.character(asignacion_taxonomica$completo), function(x){strsplit(x, ";")[[1]][2]})
asignacion_taxonomica$Class <- sapply(as.character(asignacion_taxonomica$completo), function(x){strsplit(x, ";")[[1]][3]})
asignacion_taxonomica$Order <- sapply(as.character(asignacion_taxonomica$completo), function(x){strsplit(x, ";")[[1]][4]})
asignacion_taxonomica$Family <- sapply(as.character(asignacion_taxonomica$completo), function(x){strsplit(x, ";")[[1]][5]})
asignacion_taxonomica$Genus <- sapply(as.character(asignacion_taxonomica$completo), function(x){strsplit(x, ";")[[1]][6]})
asignacion_taxonomica$Species <- sapply(as.character(asignacion_taxonomica$completo), function(x){strsplit(x, ";")[[1]][7]})
  1. Visualización de la tabla
View(asignacion_taxonomica)

Comprimir la matriz a distintos niveles taxonómicos

Usaremos la funcion taxa_levels() que cargamos previamente



taxonomic_levels("data/TablaFinal.txt", "data")
dir("data")

Ejercicio 7: Normalización de matrices de abundancia

Diapositiva 18

View(tabla)
  1. Sacar conteos totales por muestra


totales <- colSums(tabla)
totales
  1. Generar una nueva matriz












tabla_rel <- matrix(nrow=dim(tabla)[1], ncol=dim(tabla)[2])
colnames(tabla_rel) <- colnames(tabla)
rownames(tabla_rel) <- rownames(tabla)

for (i in 1:dim(tabla)[1]){
	for (j in 1:dim(tabla)[2]){
		tabla_rel[i,j] <- (tabla[i,j]/totales[j])*100
	}
}

View(tabla_rel)
colSums(tabla_rel)
  1. Escribir la matrix normalizada en nuestros archivos
exportMat(tabla_rel, file ="data/percent_integrated_matrix.txt")
  1. Correr la función abund_relativa() para los niveles de Género y Clase



abund_relativa("data/Genus.txt", "data/percent_Genus.txt")
abund_relativa("data/Class.txt", "data/percent_Class.txt")
dir("data")

Ejercicio 8: Generación de barras apiladas con Cortes de abundancia

Diapositivas: 19- 21

  1. Cargar matriz a nivel de clase en abundancias relativas



tabla <- read.delim("data/percent_Class.txt", header=T, row.names=1)
View(tabla)
dim(tabla)
  1. Determinar la media de abundancia relativa para cada Taxa (fila)



maxab <- apply(tabla, 1, mean)
length(maxab)
head(maxab)
  1. Determinar un corte de abundancia, aquí marcamos un corte al 1%
cut <- 1
  1. Obtener los taxa que tienen una abundancia mayor al corte que estamos determinando
abundancias <- maxab[which(maxab > cut)]

Número de taxas con los que nos quedaremos

length(abundancias)
  1. Generar matriz con el corte


subtable <- tabla[which(rownames(tabla) %in% names(abundancias)),]
dim(subtable)
  1. Ordenar la matrix de acuerdo a la media de abundancia relativa




tabla_taxas <- data.frame(abundancias)
tabla_taxas$taxa <- names(abundancias)
df_ordenado <- tabla_taxas[order(tabla_taxas$abundancias, decreasing =T ),]
orden <- df_ordenado$taxa
  1. Generar una nueva matriz, ya ordenada



new <- matrix(nrow=dim(subtable)[1], ncol=dim(subtable)[2])
colnames(new) <- colnames(subtable)
rownames(new) <- orden

Llenar la nueva matrix de acuerdo al orden










for (i in 1:dim(new)[2]){
	numCol <- which(colnames(tabla)==colnames(new)[i])
	for (j in 1:dim(new)[1]){
		numRow <- which(rownames(tabla)==rownames(new)[j])
		new[j,i] <- tabla[numRow, numCol]
	}
}

View(new)
  1. Normalizar al 100%. Agregar una fila con los valores de la abundancia relativa menor al corte que generamos









totales <- colSums(new)

subtable2 <- as.data.frame(matrix(ncol=dim(new)[2], nrow=1))
colnames(subtable2) <- colnames(new)
rownames(subtable2) <- paste("Abundance <", cut, "%", sep="")
for (i in 1:dim(subtable2)[2]){
	subtable2[1,i] <- 100-totales[i]
}
new <- rbind(new, subtable2)
  1. Escribir nueva matriz de abundancia
exportMat(as.matrix(new), file ="data/percent_Class_more_1_perc.txt")
  1. Contraer la información de la matrix en dos columnas (Fusionar la matriz)
melt_subtable <- melt(new)

Añadir los nombres de los taxa en una columna extra

melt_subtable$Taxa <-  as.character(rep(rownames(new), dim(new)[2]))

Ordene los nombres de los taxa de acuerdo al orden establecido




melt_subtable$Taxa <- factor(melt_subtable$Taxa, levels=rownames(new))
colnames(melt_subtable) <- c("Sample","Freq", "Taxa")
head(melt_subtable)
  1. Generar un vector de colores al azar
    Quitar la escala de grises del vector de colores
color = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]

Generar un vector de colores al azar con un color gris para los taxa menos abundantes

my_color=c(sample(color, dim(new)[1]-1, replace=FALSE), "gray")
  1. Generar mi grafica de barras apiladas con el corte de abundancia











ggplot(melt_subtable, aes(x= Sample, y=Freq, fill=Taxa)) + 
geom_bar(stat="identity") +  
	scale_fill_manual(values=my_color, limits=rownames(new)) + 
	theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1, size=6), 
	      panel.background= element_rect(fill = NA, colour = "black"), 
	      panel.grid.major = element_line(colour = "grey90"), 
	      legend.key = element_rect(fill = "white"), 
	      legend.text=element_text(size=6), 
	      legend.position="bottom") + 
	labs(x = "", y = "Relative abundance (%)", fill="") + 
    scale_y_continuous(breaks=c(0, 20, 40, 60, 80, 100))

Guardar imagen con función ggsave()



ggsave("figuras/sock_plots_Class_more_1_perc.png", 
       device = "png", dpi = 300, units = "px", bg = "white")

Ejercicio 9: Library qiime2R

Leyendo y manejando artifacs en R

Diapositivas: 23-29

  1. Cargar librerias de trabajo


library(qiime2R)
library(tidyverse)
  1. Leer el "artifac" que genero QIIME2 (table.qza)
tabla_qza <- read_qza("qiime2/table.qza")
  1. Visualizar la estructura de nuestro objeto


str(tabla_qza)
class(tabla_qza)
  1. Extraer la información del objeto tabla_qza como una lista
names(tabla_qza)
  1. Acceder a la información de las features (ASVs)


View(tabla_qza$data)
head(tabla_qza$data)
  1. Asociar la taxonomía



taxonomia_qza <- read_qza("qiime2/taxonomy_blast.qza")
str(taxonomia_qza)
class(taxonomia_qza)
  1. Acceder a la taxonomía
View(taxonomia_qza$data)
  1. Dividir en niveles taxonómicos la taxonomía usando parse_taxonomy()


taxa_sep <- parse_taxonomy(taxonomia_qza$data)
View(taxa_sep)
  1. Comprimir el data.frame a un nivel taxonómico en particular (Nivel de clase)


class_table <- summarize_taxa(tabla_qza$data, taxa_sep)$Class
View(class_table)

Gráfica de barras apilada y heatmaps usando la información de metadata

  1. Cargar el archivo de metadata.tsv


metadata <- read_q2metadata("qiime2/metadata.tsv")
View(metadata)
  1. Generación de gráfica de barras apiladas
taxa_barplot(class_table, metadata,"Zone")
  1. Generación de gráfica de Heatmap
taxa_heatmap(class_table, metadata,"Type")
  1. Cambiar gamma de colores de acuerdo a la gramática de ggplot2











taxa_barplot(class_table, metadata,"Zone") +
  scale_fill_manual(values=my_color) +
  theme(text = element_text(size=10), legend.position="bottom") +
  scale_y_continuous(breaks=c(0, 20, 40, 60, 80, 100)) +
  labs(fill="")

taxa_heatmap(class_table, metadata,"Type") + 
  scale_fill_gradientn(colours = terrain.colors(10))

taxa_heatmap(class_table, metadata,"Type") +
  scale_fill_gradient(low = "purple", high = "green3", na.value = NA)
  1. Guardar imagen en png



png("figuras/taxa_heatmap.png", width = 10*300, height = 8*300, res = 300, pointsize = 8)
taxa_heatmap(class_table, metadata,"Type")
dev.off()

Generación de un archivo de tipo phyloseq










physeq_qiime2 <- qza_to_phyloseq (
  features="qiime2/table.qza",
  tree="qiime2/rooted-tree.qza",
  taxonomy="qiime2/taxonomy_blast.qza",
  metadata = "qiime2/metadata.tsv")

physeq_qiime2
class(physeq_qiime2)
str(physeq_qiime2)

Alfa y Beta diversidad

Ejercicio 10: Cálculo de los estimadores de alfa diversidad usando phyloseq

Diapositivas: 30-34

Cargar los datos











variables <- read.delim("data/abiotic_variables.txt", header=T, row.names=1)
abun_table <- read.delim("data/Genus.txt", header=T, row.names=1) 
abun_table_mat <- as.matrix(abun_table)
dim(abun_table)
View(abun_table)

taxrank_mat <- as.matrix(rownames(abun_table_mat))
colnames(taxrank_mat) <- ("TAX")
rownames(taxrank_mat) <- rownames(abun_table_mat)
head(taxrank_mat)

Generar objetos de tipo phyloseq con las funciones de phyloseq especificas

  1. Matriz de OTUs (En nuestro caso es la matriz de Géneros)
OTU <- otu_table(abun_table_mat, taxa_are_rows=TRUE)
  1. Una matriz con las anotaciones de los OTUS (En nuestro caso es la anotación taxonómica a nivel de género)
TAX <- tax_table(taxrank_mat)
  1. Cargar tabla de metadatos
MET <- sample_data(variables)
  1. Construimos nuestro objeto de phyloseq


physeq <- phyloseq(OTU, TAX, MET)
str(physeq)
  1. Obtenemos los índices de diversidad y los colocamos en una tabla


diversity_index <- data.frame(round(estimate_richness(physeq, measures=c("Observed", "Chao1", "Shannon", "Simpson")), 2))
diversity_index

Usar el objeto que construimos con qiime2R



diversity_index2 <- data.frame(round(estimate_richness(physeq_qiime2, measures=c("Observed", "Chao1", "Shannon", "Simpson")), 2))
diversity_index2
  1. Escribimos la tabla
write.table(diversity_index, "data/alpha_diversity_index.txt", quote = FALSE, sep = "\t", row.names = TRUE, col.names=TRUE)
  1. Generar un boxplot usando el índice de Shannon










diversity_index$Clase <- sapply(as.character(rownames(diversity_index)), 
                                function(x){strsplit(x, "_")[[1]][1]})

ggplot(diversity_index, aes(x=Clase, y=Shannon, fill=Clase)) +
  geom_boxplot(show.legend=FALSE) + 
  geom_point(color="black", alpha=0.9, position = "jitter", shape=21, show.legend=FALSE) +
  labs(x = "", y = "Indice de Shannon", fill="") +
  scale_fill_manual(values=c("royalblue2", "tomato2", "seagreen3")) +
  theme_minimal()
ggsave("figuras/boxplot_shannon.png", device = "png", dpi = 300, units = "px", bg = "white")

Ejercicio 11: Generación de PCoA/NMDS con phyloseq

Diapositivas: 35-38

Generación del PCoA con el método de Bray-Curtis usando la función ordinate()



ordination <- ordinate(physeq, "PCoA", "bray")
plot_ordination(physeq, ordination, color="grupo")

Podemos usar la información de alguna variable continua

plot_ordination(physeq, ordination, color="Deep")

Cambiarle el color a otra paleta de colores



plot_ordination(physeq, ordination, color="Deep") + 
  scale_colour_gradientn(colours= topo.colors(10))

Vamos a usar la sintáxis de ggplot2 para mejorar las gráficas






plot_ordination(physeq, ordination, color="grupo") + 
  stat_ellipse() +
  scale_color_manual(values= c("deeppink1", "darkturquoise", "gold")) +
  labs(title="PCoA") + 
  theme_bw()

Ejemplo de como generar un PCoA cuando tenemos la información del árbol filogenético

plot_ordination(physeq_qiime2, ordinate(physeq_qiime2, "PCoA", "wunifrac"), color="Type")

Guardar nuestro gráfico





plot_ordination(physeq, ordination, color="grupo") + stat_ellipse() +
  scale_color_manual(values= c("deeppink1", "darkturquoise", "gold")) +
  labs(title="PCoA") + theme_bw()
ggsave("figuras/PCoA_bray_points_phyloseq.png", device = "png", dpi = 300, units = "px")

Generación de un NMDS





plot_ordination(physeq, ordinate(physeq, "NMDS", "bray"), color="grupo") +
  scale_color_manual(values= c("deeppink1", "darkturquoise", "gold")) +
  labs(title="NMDS") + theme_bw()
ggsave("figuras/NMDS_bray_points_phyloseq.png", device = "png", dpi = 300, units = "px")

Ejercicio 12: ANOSIM y ADONIS

Diapositivas: 38-39

Normalizar usando el método de metagenomeSeq.

  1. Cargar matriz de conteos absolutos
datos <- loadMeta("data/Genus.txt")
  1. Transformarla a un objeto de metagenomeSeq llamado "MRexperiment"
MRexp <- newMRexperiment(datos$counts)
  1. Calcular el percentil para el cual sumar los conteos absolutos
p <- cumNormStatFast(MRexp)
  1. Normalizamos la matriz


normalized_matrix <- cumNormMat(MRexp, p = p)
head(normalized_matrix)
  1. Escribimos la matriz normalizada
exportMat(normalized_matrix, file ="data/normal_Genus.txt")
  1. Transponer la matriz normalizada para que las filas sean las muestras y los taxa las columnas
transp_mat <- t(normalized_matrix)
  1. Generar grupos




muestras <- data.frame(estacion=rownames(transp_mat))
muestras$clase <- sapply(as.character(muestras$estacion), 
                         function(x){strsplit(x, "_")[[1]][1]})
grupos <- muestras$clase

ANOSIM (Diferencias entre todos los grupos)







str(transp_mat)
resultado_anosim <- anosim(transp_mat, 
                           grouping= grupos, 
                           permutations = 999, 
                           distance = "bray")
resultado_anosim

ADONIS pareado (comparación pareada de los grupos)





pairwise.adonis(x= transp_mat, grupos, 
                sim.function = 'vegdist', 
                sim.method='bray', 
                p.adjust.m='bonferroni')

Información de la sesión y Listado de objetos



ls()
sessionInfo()

Información adicional y enlaces

Gamma de Colores en R
Colores HTML
Paleta de colores Opción 1
Paleta de colores Opción 2
Paleta de colores Opción 3
metagenomeSeq
phyloseq
qiime2R

tags: UUASMB Bitácoras Metagenómica
Last changed by 
 
Tina Godoy
Head of the Department of Bioinformatics in Infectious Diseases at the National Institute of Public Health in Mexico.
0
270

Read more

Metagenómica funcional 2024. Tema 4B

Centro de Investigaciones Biológicas del Noroeste S.C. La Paz, Baja California Sur.Noviembre 2024

Nov 21, 2024
Metagenómica funcional 2024. Tema 2

downloadsbuildchat

Nov 19, 2024
Metagenómica funcional 2024. Tema 3

Centro de Investigaciones Biológicas del Noroeste S.C. La Paz, Baja California Sur.Noviembre 2024

Nov 19, 2024
Genómica Comparativa

downloadsbuildchat

Nov 13, 2024
Read more from Tina Godoy
Published on HackMD