ampvis2

Un buen programa para realizar varias tipos de visualizaciones es ampvis2 que trabaja en R. Necesitamos dos archivos, una tabla de OTUs y una de metadatos y tener instalado R y de preferencia RStudio.

mg_classifier a partir de la versión 1.8.0 ya genera una tabla de otus compatible (OTU_table.tsv) y la tabla de metadatos hay que crearla según las indicaciones aquí mencionadas. Es importante que los códigos de las muestras en ambas tablas sean idénticos.

Instalación de ampvis2

En la consola de R, dentro de RStudio pegar los siguientes comandos:

install.packages("remotes")
remotes::install_github("MadsAlbertsen/ampvis2")

Formato de las tablas

Tabla de OTUs

La tabla de OTUs tiene que tener ciertas características, ser una tabla delimitada por tabuladores, pero principalmente al final tener las columnas de la taxonomía de cada OTU con un nombre concreto y solo los 7 niveles clásicos en inglés, empezando por Kingdom.

Ejemplo

OTU F3D1 F3D141 F3D142 Kingdom Phylum Class Order Family Genus SpeciesOTU_3 2 0 0 Bacteria Bacteroidetes Bacteroidia Bacteroidales S24-7_f Unclassified_g Unclassified_sOTU_5 0 0 0 Bacteria Firmicutes Clostridia Clostridiales Lachnospiraceae Unclassified_g Unclassified_sOTU_6 0 0 0 Bacteria Bacteroidetes Bacteroidia Bacteroidales S24-7_f Unclassified_g Unclassified_sOTU_7 0 0 0 Bacteria Firmicutes Clostridia Clostridiales Ruminococcaceae Unclassified_g Unclassified_sOTU_8 0 0 0 Bacteria Firmicutes Clostridia Clostridiales Ruminococcaceae GQ175418_g Unclassified_s

Tabla de metadatos

La tabla de metadatos tiene un formato semejante, igual delimitada por tabuladores:

sample  dpw time
F3D1    1   Early
F3D141  141 Late
F3D142  142 Late

El único requisito es que el nombre de la muestra (sample) debe ser igual al que tiene la tabla de OTUs.

Los datos aquí usados se pueden bajar de las siguientes ligas:

Cargar los archivos a RStudio

Importar tabla de OTUs

Dar click en Import Dataset y luego en "From Text..."

Seleccionar el archivo de tabla de otus, en este ejemplo, OTU_table.tsv

Revisar que este marcado Heading: Yes y Row names: Use first column

Ahora ya aparece la tabla de otus cargada en el panel superior derecho (flecha roja).

Importación de metadatos

Volvemos a hacer lo mismo para importar los metadatos.

Cargar el paquete de ampvis2

Abrir la pestaña Packages y dar click en el cuadro de ampvis2

Crear una objeto en ampvis2

Ya en la consola hay que crear un matriz (aquí llamada SOP) con los datos cargados

> SOP <- amp_load(otutable = OTU_table, metadata = metadata)

Se creará una nueva matriz con el nombre SOP (Fig Der), ésta matriz es la que usaremos para generar los análisis y visualizaciones. Como se ve en el ejemplo, al crear el objeto lanza una advertencia que hay una inconsistencia entre la tabla de OTUs y la de metadatos, esto es por que no hay la muestra llamada "Mock" en la tabla de metadatos, por lo que no utilizará esos datos que sí están en la tabla de OTUs.

Análisis

Rarefacción

Quizá el primer análisis a realizar es una curva de rarefacción con la función amp_rarecurve:

> amp_rarecurve(data = SOP, stepsize = 50, color_by = "time")

Como se puede apreciar en la figura de la abajo a la izquierda, se generó una curva de rarefacción coloreada por la columna "time" de los metadatos. El comando stepsize implica que tan frecuente hace el cálculo para generar la gráfica, a menor número, mas fina es la curva. Podemos mejorar la gráfica agrupando las muestras por algún elemento de los metadatos, en este caso de nuevo por "time".

> amp_rarecurve(data = SOP, stepsize = 50, color_by = "time", facet_by = "time")

El resultado está en la figura de abajo a la derecha.

Boxplot

Podemos visualizar los datos como un boxplot:

> amp_boxplot(data = SOP, group_by = "time", tax_add = "Phylum")

Al igual que en el caso anterior, están agrupadas por la columna "time" de los metadatos y a nivel género, que es el default, pero se puede cambiar. A cada género se le añadió además el phylum al que pertenece con tax_add.

Análisis de Componentes Principales (PCA)

Para obtener un PCA:

> amp_ordinate(data = SOP, sample_color_by = "time")

pero también se le puede agregar un marcho a los puntos (Fig. abajo Der.):

> amp_ordinate(data = SOP, sample_color_by = "time", sample_colorframe = "time")

Heatmap

Para obtener un heatmap agrupados por "time" (group_by = "time"), a nivel género (tax_aggregate) pero mostrando también el phylum (tax_add = "Phylum") y mostrando los 20 taxa mas abundantes (tax_show = 20):

> amp_heatmap(data = SOP, group_by = "time", tax_show = 20, tax_aggregate = "Genus", tax_add = "Phylum")

Diagram de Venn

Podemos crear un diagrama de Venn basándonos en un parámetro de los metadatos (group_by = "time") y a cierto porcentaje considerado como las OTUs presentes en el núcleo o core (cut_f = 95), en este caso el 95%:

> amp_venn(data = SOP, group_by = "time", cut_f = 95)


Formateado de las gráficas

Debido a que estas gráficas son creadas con ggplot, es posible embellecerlas más. Por ejemplo, si usamos la función para la creación de una curva de rarefacción simple:

amp_rarecurve(data = SOP, stepsize = 50, color_by = "STAGE", facet_by = "STAGE")

se generará la figura de abajo a la izquierda, pero si agregamos los comandos de ggplot ( + theme()) crearemos una mejor figura (abajo derecha), las letras son mas grandes, sin leyenda, líneas mas gruesas, etc.

amp_rarecurve(data = SOP, stepsize = 50, color_by = "time", facet_by = "time") + 
    theme(
        axis.title.x = element_text(size = rel(1.4)),
        axis.text.x = element_text(size = 14, angle = 0),
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(size = rel(1.4)),
        axis.text.y = element_text(size = 14),
        axis.ticks.y = element_blank(),
        legend.position = "none",
        strip.text = element_text(size = 14, face = "bold"),
        panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"),
    )  + geom_line(size = 1.0)

Lo mismo podemos hacer para la gráfica de componentes principales (PCA), el código simple (abajo Izq):

amp_ordinate(data = SOP, sample_color_by = "time", sample_point_size = 4, sample_label_by = "sample", sample_label_size = 5)

Aún mejor incluyendo los principales géneros (Abajo Der):

amp_ordinate(data = SOP, sample_color_by = "time", sample_point_size = 4, sample_label_by = "sample", sample_label_size = 5, species_plot = TRUE, species_nlabels = 16, filter_species = 1, species_label_size = 4, sample_colorframe = "time")

y el código tuneado (hasta abajo Der):

amp_ordinate(data = SOP, sample_color_by = "time", sample_point_size = 4, sample_label_by = "sample", sample_label_size = 5, species_plot = TRUE, species_nlabels = 16, filter_species = 1, species_label_size = 4, sample_colorframe = "time") +
    theme(
        axis.title.x = element_text(size = rel(1.2)),
        axis.title.y = element_text(size = rel(1.2)),
        axis.text.x = element_text(size = 14, vjust = 1), 
        axis.text.y = element_text(size = 14),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        strip.text = element_text(size = 14), # facets
        panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14, face = "bold"),
    )

Heatmap mejorado y con escalas de grises

> amp_heatmap(data = SOP, group_by = "time", facet_by = "time", tax_aggregate = "Genus", tax_add = "Phylum", tax_show = 20, round = 2, plot_values_size = 6, color_vector = c("#f2f2f2", "#909190")) +
  theme(
  axis.text.x = element_text(size = 14, angle = 45, vjust = 1),
  axis.text.y = element_text(size = 14, face = "italic"),
  axis.ticks.x = element_blank(),
  axis.ticks.y = element_blank(),
  strip.text = element_text(size = 14, face = "bold"),
 )