# Datos 16S restricción calórica - round 2
###### tags: `16S rRNA`
[ToC]

## Análisis exploratorio inicial
:::info
Estos gráficos son **PCAs** hechas con la configuración por defecto de Phyloseq, que usa la **distancia Bray-Curtis**.
:::

:::danger

:::

Incluyendo más muestras, nos sale esto (el panel de arriba colorea por grupos, y el panel de abajo solo distingue los grupos con restricción calórica o alimentación _ad libitum_ _**a tiempo 0**_, es decir, el G7 frente al resto):


Aquí incluyo un par de PCs más, aunque ya no tienen mucha información:

Vemos que sí hay grupos que se separan, pero ojo porque aquí mezclamos varios tiempos distintos. Quedándonos solo con las muestras a tiempo 0:


> En el desplegable dejo una PCA separando las muestras por tiempos:
> :::spoiler
> 
> :::

Seguimos viendo que los puntos tienen cierta agrupación, sobre todo los grupos 1 y 2 en los paneles con PC1/PC3 y PC2/PC3. Los ratones en RC (que a tiempo 0 es solo el grupo 7) sí parecen quedar apartados cuando analizamos los dos primeros componentes principales.

### Grupos control (1 y 2)
Por otro lado, dijimos que nos interesaría ver los ratones de los grupos 1 y 2 en los tres tiempos, esperando verlos juntos en tiempo 0, algo separados en tiempo 1, y por ver qué pasa en tiempo 2:


Y esto es justo lo que vemos :-) solo sale significativo el t1, con Pr(>F) = 0.001 (estadística: PERMANOVA con 999 permutaciones)

:::info
Comentando esto en una reunión que tuvimos, parece interesante que en tiempo 2 se nos mezclen los ratones de ambos grupos. Durante mi TFM dimos bastantes vueltas al llamado "principio de Anna Karenina". Este principio se llama así por la frase de _"todas las familias felices se parecen unas a otras, pero cada familia infeliz lo es a su manera"_ del libro de Anna Karenina, que aplicado a microbioma sería algo como que todos los microbiomas sanos se parecen, pero cada microbioma "no sano" lo es a su forma...
Lo que sugiere este principio es que el paso de un microbioma sano a uno no sano se basa sobre todo en cambios estocásticos, más que deterministas, de forma que no se podría predecir cómo una enfermedad (en este caso, el cáncer) va a afectar al microbioma de un individuo (o ratón), porque en cada uno va a ser distinta. Esto se ha visto en ecosistemas y en algunos estudios que miran cambios en la microbiota en enfermedades. Viendo estos resultados, pensamos que podíamos dejar esto aquí por si os interesaba. Un par de artículos sobre ello:
- Zaneveld, J., McMinds, R. & Vega Thurber, R. Stress and stability: applying the Anna Karenina principle to animal microbiomes. Nat Microbiol 2, 17121 (2017). https://doi.org/10.1038/nmicrobiol.2017.121
- Ma ZS. Testing the Anna Karenina Principle in Human Microbiome-Associated Diseases. IScience 2020;23:101007. https://doi.org/10.1016/j.isci.2020.101007
:::

:::warning
Probablemente queramos hacer más PCAs a lo largo de los análisis, pero primero vamos a seguir y ya iremos dibujando a medida que nos encontremos con las preguntas :-)
:::

## Análisis de alfa-diversidad
:::info
Aquí había que solucionar la forma de hacer las gráficas, para que sean más interpretables, y añadir los grupos 6 y 7 para poder usarlos como referencia.
:::

### Análisis longitudinal
Habíamos quedado en imitar las figuras que creé para _Akkermansia_ y _Lactobacillus_, con lo que nos quedaría una gráfica así:


En este caso, no salen cambios significativos a lo largo del tiempo en ninguno de los casos (estadística: Wilcoxon pareado + corrección FDR). De hecho, si hiciéramos el mismo gráfico pero en lugar de con boxplots, uniendo con líneas los puntos del mismo ratón, veremos que cada uno hace una cosa diferente:


He calculado las varianzas de cada grupo en cada punto temporal. Aquí podemos comprobar que los grupos 6 y 7 tienen una varianza mucho mayor que los demás, como ya intuíamos por las gráficas:


En esta gráfica tenemos la evolución de la alfa-diversidad a nivel individual, pero habíamos dicho que además podría estar bien mostrar en un mismo plot cómo evoluciona la media de alfa-diversidad en cada grupo:


Mostrando las barras de error (intervalo de confianza 95%):


Estas son las visualizaciones que considero que han quedado mejor. Debajo pongo un desplegable con algunas pruebas que he hecho, para que queden aquí y por si consideramos que puede ser interesante avanzar con alguna de ellas. Esto incluye las mismas gráficas poniendo además los puntos individuales de cada ratón, que era la idea inicial pero no conseguía hacer que se vieran bien.

:::spoiler
Mostrando los puntos además de las medias:

Otra posibilidad es colorear el área que dejarían las barras de error:
Pudiendo poner los puntos por encima:

Y las mismas gráficas, pero incluyendo además los grupos 6 y 7:

:::

:::danger
Lo que me está dando miedo es que tengamos tanta variabilidad en los ratones de los grupos 6 y 7 que no nos sirvan de mucho para las comparaciones
:::

### Análisis en cada tiempo
Hasta aquí tendríamos el análisis **longitudinal**, pero nos sigue quedando la comparación entre grupos para cada uno de los tiempos. Esta gráfica la tenemos, solo nos faltaba añadir los grupos 6 y 7 (donantes):


:::danger
**Duda:** Aquí no sé si la corrección de p-values debería hacerse dentro de cada medida (como está hecho aquí, se calculan los p-values de Chao1 y se ajustan, luego los de Shannon, y luego los de Simpson), o al final (se calculan los p-values de las tres medidas y se ajustan todos de golpe). Antes lo tenía hecho de esta segunda forma y por eso salían menos p-valores significativos.
:::

## Análisis de beta-diversidad
### Restricción calórica _vs._ alimentación _ad libitum_
:::info
Recordemos que en realidad esto es muy parecido a lo que hacíamos al principio con las PCAs. En ambos casos estamos calculando alguna medida de distancia entre muestras y después usando gráficas de ordenación para visualizarlo. Antes usamos distancia Bray-Curtis + PCA, aquí usamos distancia Weighted UniFrac + MDS.
:::

Aquí queríamos también ver los grupos 1 y 2:


>Haciendo el adonis sale diferencia solo en tiempo 1, con Pr(>F) = 0.001 (estadística: PERMANOVA con 999 permutaciones).

Si añadimos los grupos 6 y 7 y los comparamos con el t1 de los grupos 1 y 2 (cuando podríamos ver el efecto la dieta), esperaríamos ver que el 1 y el 6 se parecen (AL) y lo mismo con el 2 y el 7 (RC), pero lo que nos encontramos es que los ratones del G7 están muy alejados del resto y desplazan el resto de la gráfica:


También queríamos ver, en todos los tiempos, si hay diferencias entre los animales en restricción calórica / con trasplante de restricción calórica y los de *ad libitum*.

## Análisis de dispersión
### Dispersión de todos los ratones en cada tiempo
Recordemos esta gráfica que habíamos visto, en la que mostramos la beta-diversidad separando los ratones por puntos:


Nos había llamado la atención que parece que en tiempo 1 las muestras son menos variables que en los otros dos tiempos. Vamos a analizar si esto es así. Para ello utilizaremos la función `betadisper` del paquete `vegan`:


Esta figura muestra todos los ratones y su distancia al centroide, que es la mediana de los puntos en el espacio (más sobre esto en la ayuda de la función). Los centroides se han calculado agrupando por tiempo: tenemos uno para todas las muestras a tiempo 0, otro para tiempo 1 y otro para tiempo 2. Para comprobar si realmente hay diferencias en las **distancias de los puntos al centroide de cada grupo**, es decir, si hay diferencia en la dispersión de los puntos, lo miramos con `permutest`, dándole el parámetro `pairwise = TRUE` para que haga también las comparaciones dos a dos y nos diga qué grupos tienen diferencias significativas (si es el caso):

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups      2 0.03746 0.0187322 4.4957    999  0.012 *
Residuals 166 0.69167 0.0041667                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          t0        t1    t2
t0           0.1030000 0.005
t1 0.1070789           0.171
t2 0.0041112 0.1514030   

Vemos que nos sale un p-valor significativo para los grupos t0 y t2, pero no para ninguna comparación con t1.

:::info
De aquí lo que sacamos es que los puntos en tiempo final tienen significativamente ___menor dispersión___ que los puntos en tiempo inicial.

> He sacado la misma figura con ggplot porque quería modificarla un poco (y así además tenemos más control sobre el plot si queremos cambiar cosas a futuro):
> :::spoiler
> 
> :::
:::

Me surgió la duda de si los grupos 6 y 7, que hemos visto por ejemplo en la alfa diversidad que eran muy variables (más que el resto), podrían estar afectando a este resultado, y que realmente con los ratones de los grupos 1-5, que son los que llegan hasta tiempo final, no veamos esa diferencia. Probé a correr este mismo análisis excluyendo los ratones de los dos grupos donantes y me salió lo mismmo (es decir, diferencias significativas entre t0 y t2). El p-valor es algo más alto, pero sigue siendo < 0.05 (da 0.045).

A raíz de esto pensé: ¿y si el grupo 6 y 7 realmente tienen mayor dispersión que el resto? Esta es la pregunta que intento resolver en el siguiente apartado.

### Dispersión de los ratones de cada grupo
Ahora vamos a mirar si vemos diferencias en la dispersión de los puntos de los ratoncitos pertenecientes a distintos grupos (estratificando de nuevo por tiempos). Así también podemos ver si esa variabilidad que veíamos en las medidas de alfa-diversidad se da también cuando miramos beta-diversidad.

En este caso pongo directamente las figuras de `ggplot`, porque creo que son las que mejor se ven. Los puntos con borde negro son los centroides, y los más pequeños sin borde son las muestras:


> Las figuras que saca `betadisper` las dejo aquí debajo (son iguales, pero con otros colores y formas):
> :::spoiler
> **tiempo 0:**
> 
> **tiempo 1:**
> 
> **tiempo 2:**
> 
> :::

:::info
En ninguno de los tiempos salen diferencias significativas en la dispersión de los puntos entre distintos grupos. **Recordemos que la pregunta que nos estamos haciendo aquí es esta, sobre la _dispersión_, no si al representar los puntos basados en la distancia Weighted UniFrac vemos diferencias.**
:::

## Por ratoncitos
Dejo este apartado porque originalmente todo el análisis con `betadisper` venía motivado de haber comentado estas figuras en la reunión:

:::spoiler
:::

Las dos figuras representan lo mismo (beta-diversidad uniendo los puntos de un mismo ratón en un triángulo), solo que en una mostramos todo a la vez y en otra separamos por grupos. Lo que habíamos comentado es que podría ser interesante ver si los triángulos de los ratones de un grupo tienden a ser más grandes de los de otro, cuando miramos todos los puntos a la vez.

> En concreto, habíamos dicho que esperaríamos que los triágulos de los ratones en RC fueran más grandes que los de los ratones AL; porque entre tiempo 0-1 esperamos que la microbiota de los ratones RC cambie más que las de los ratones AL (si lo tengo bien apuntado).

La cosa es que al hacer el análisis de `betadisper` separando por grupos no estamos viendo exactamente eso: estamos viendo variabilidad _entre_ ratones, no _intra_ ratón, que era la pregunta que teníamos originalmente.

Quizá una buena aproximación sea ver todos los puntos y seaprar por grupos (es decir juntar los tres paneles de la figura anterior), aunque no sé si es exactamente lo que queremos:

> Aquí en cada grupo estamos mezclando muestras de t0, t1 y t2. No hay diferencias significativas en la dispersión de los datos de cada grupo.

:::danger
Esto es parecido a lo que yo creo que habría que hacer:
> * Para ver esto, en `betadisper` tenemos que usar el ID de ratón, no el grupo de dieta ni el tiempo. El tema es que al hacer esto, nos compara todos los ratones entre sí: en lugar de mirar si en general los ratones del grupo 1 tienen más dispersión que los del 2, veremos si el ratón 153 tiene más dispersión que el 212 (no sale nada significativo).
> * Creo que juntar todos los datos de t0, t1 y t2 y ver el área por grupos no es exactamente lo mismo, porque lo que queríamos es ver a nivel de cada ratón individualmente (no sé si le estoy dando demasiadas vueltas a esto, probablemente sí)
> * Lo mejor que he conseguido es mirar en cada ratón cuáles son las distancias de los puntos al centroide... pero esto sigue sin ser el área del triángulo que se forma:
> 

Esto no me convence nada y creo que me estoy complicando la vida.
:::

## Géneros de interés
### *Akkermansia* y *Lactobacillus*
:::info
Aquí quedaba:
- [x] Añadir los grupos 6 y 7 a las gráficas
- [x] Aprovechando esto, queremos **comprobar si estos grupos tienen *Lactobacillus***
- [ ] Repetir las gráficas mostrando _fold change_
:::
#### *Akkermansia*

### *Lactobacillus*

> Extra (para ver distribución de los puntos, creo que puede venir bien porque en las anteriores figuras no nos hacemos a la idea de cuántos puntos hay en el 0):
> :::spoiler
> 

### Otros géneros
:::info
- [ ] Aquí habíamos dicho de **repetir el ANCOM** con las preguntas que nos interesan.
:::

## Más cosas pendientes
:::info
- [ ] Curvas densidad grupo tPBSr
:::

Esta figura la habíamos planteado para ver si había mucho desplazamiento entre la curva de tiempo 0 y la de tiempo 1, para ver cuántos géneros se recuperan tras el FMT por sí solos (ya que aquí no hay ningún aporte de microbiota, solo PBS). El tema es que tenemos tantos géneros con abundancias próximas a 0 que no queda una figura muy ilustrativa:

> He puesto las barritas una junto a otra en lugar de superpuestas porque si no, como tienen alturas tan parecidas, no se ve nada.

:::spoiler
La misma figura pero con la línea de densidad en lugar del histograma:
:::

:::info
Quizá otra aproximación pueda ser hacer un heatmap en el que mostremos cómo cambian las especies entre tiempo 0 y tiempo 1 (tipo una matriz de correlación, no sé si como tal o mejor mostrando simplemente fold change, o algo así?)
:::

:::info
- [ ] Perfiles k-means
:::