# ML y Microbioma ###### tags: `Estudiantes` `ML` ### Informe (borrador) * [Link](https://docs.google.com/document/d/1tkgc5EabgwfVQuJzDWC2Cd7omFBV2FdiGKwQeAs3tcE/edit?usp=sharing) del informe sobre Machine Learning. * [Link](https://docs.google.com/document/d/1RsZGjHIWmo8GaJx_r2jX_fViL8n1lNQU8xncEG8n5K8/edit?usp=sharing) del software de ML traducido al inglés. * [Link](https://docs.google.com/spreadsheets/d/18rC8WMnOdUJ89SYRukASokPCRUve7XDZzq2FBSxAQhU/edit?usp=sharing) del Excel donde iré añadiendo información sobre software de ML y microbiota que vaya encontrando. * [Tabla](https://docs.google.com/spreadsheets/d/1BoDhS1nNIM5p2RrR3G6u-x52wy1no7ev5qUGEh-RpJ4/edit?usp=sharing) con artículos para revisión de métodos de pre-procesado y transformaciones (revisión COST Action) ### Papers * [Revisión](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7962872/) de la COST Action. ### Picture Your Microbes * [Tabla](https://docs.google.com/spreadsheets/d/1p14nwGE2Q9mUAXRpzaqOXLWY9Xd-a552/edit?usp=sharing&ouid=101078042900002591640&rtpof=true&sd=true) de metadatos de PYM. ### Trabajo con datos en qiime2 Tras iniciar sesión en Gandalf, creé un nuevo conda environment: ``` conda create --name victor ``` Tras esto, en dicho environment, instalé la última versión de qiime: ``` conda activate victor wget https://data.qiime2.org/distro/core/qiime2-2022.2-py38-linux-conda.yml conda env create -n qiime2-2022.2 --file qiime2-2022.2-py38-linux-conda.yml #Optional clean-up rm qiime2-2022.2-py38-linux-conda.yml ``` En la ruta ~/victor/PYM, existen dos zips correspondientes a los datos: PYM-1.zip y PYM-2.zip. En primer lugar, vamos a unificar los datos en un único archivo, para lo cual es necesario descomprimir ambos ficheros .zip: ``` unzip PYM-1.zip unzip PYM-2.zip ``` Una vez han sido descomprimidos ambos ficheros .zip, vamos a crear un fichero .csv, denominado "Manifest", con la localización de cada fichero y el ID de cada muestra, antes de importarlos en qiime2, y éste se llamará data.csv. Para ello he creado un script llamado "script.pl" (basado en el HackMD de Anna) con el siguiente código (info del comando vi [aquí](https://www.howtoforge.com/faq/how-to-edit-files-on-the-command-line/)). ``` touch script.pl #para crear el fichero vi script.pl #para abrir el fichero en la terminal y poder editarlo #copio y pego lo siguiente en dicho fichero con la herramienta vi: #!/usr/bin/perl use Cwd 'abs_path'; open(OUT,">data.csv"); foreach $ar (@ARGV) { #Para sacar el $name del archivo: $name=$ar; @field = split (/_/, $ar); $name2 = @field[0]; #Para sacar el $path: $path =abs_path($ar); if ($name =~ /1.fq.gz/) { print OUT "$name2,$path,forward\n"; }else{ print OUT "$name2,$path,reverse\n"; } } close OUT; ``` Creado el script, lo ejecuté: ``` perl script.pl *.fq.gz ``` Ahora, vamos a añadir una fila inicial al fichero data.csv ("Manifest") con la siguiente información: > sample-id,absolute-filepath,direction Para ello, usamos el siguiente comando: ``` vi data.csv ``` Dentro del programa vi, usando el comando i, podemos editar el archivo, insertar la línea mencionada anteriormente y luego guardar el fichero modificado pulsado ESC y luego escribiendo :x. Ya tenemos el fichero .csv con la línea inicial añadida. Tras esto, ya podemos importar los archivos en qiime2. Creo, por otro lado, que no es necesario hacer demultiplexing (para conocer qué secuencias proceden de qué muestras usando los barcodes (sirven como identificadores)), porque parece ser necesario disponer de un fichero por separado con la secuencia de los barcodes de las secuencias. Además, según [esta web](https://astrobiomike.github.io/amplicon/demultiplexing#:~:text=Demultiplexing%20refers%20to%20the%20step,had%20all%20be%20sequenced%20together.), si las muestras están separadas en archivos .fastq individuales (que es nuestro caso), no es necesario hacer demultiplexing. * **IMPORTANCIÓN DE LOS DATOS EN QIIME2** Dicho esto, importamos los datos en qiime y, con *qiime demux summarize*, obtenemos un informe de la distribución de las calidades de las secuencias en cada posición: ``` conda activate qiime2-2022.2 qiime tools import \ --type 'SampleData[PairedEndSequencesWithQuality]' \ --input-path data.csv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33 \ qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv ``` Tras ejecutar esto (basado en el HackMD de Anna), obtenemos un fichero llamado "demux.qzv", el cual podemos visualizar. Para ello, podríamos usar el siguiente comando: ``` qiime tools view demux.qzv ``` Pero no es posible al estar en Gandalf, por lo que me he pasado el fichero desde Gandalf hasta el PC local: ``` scp lmarcos@172.16.202.45:~/victor/PYM/demux.qzv . ``` Tras esto, he subido el fichero [aquí](https://view.qiime2.org/), obteniendo los siguientes resultados: ![](https://i.imgur.com/AFt838c.png) Comparando con un ejemplo que da la página web de qiime ([este](https://view.qiime2.org/visualization/?type=html&src=https%3A%2F%2Fdocs.qiime2.org%2F2022.2%2Fdata%2Ftutorials%2Fmoving-pictures%2Fdemux.qzv)), esto apunta a que los primers y las lecturas de baja calidad ya han sido eliminados. Se supone que con estos informes de calidad podría determinar los valores de los parámetros *p-trunc-len-f* y *--p-trunc-len-r* del siguiente paso (en el que se usa DADA2), pero esto parece evidenciar que los valores de estos parámetros son 0. * **PROCESADO DE LOS DATOS** Laura me confirmó posteriormente que efectivamente los primers ya habían sido eliminados, así que he ejecutado DADA2 (que sirve para llevar a cabo un control de calidad: eliminar primers y las regiones de baja calidad de las secuencias) de la forma siguiente: ``` qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-trim-left-f 0 \ --p-trim-left-r 0 \ --p-trunc-len-f 0 \ --p-trunc-len-r 0 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qza ``` Puse 0 en los parámetros porque no podía omitirlos; qiime me obliga a indicarlos. Aunque ya hayan sido eliminados los primers y las lecturas de baja calidad, tengo que hacer el paso de DADA2 porque me genera un fichero llamado *table.qza* que es necesario para el siguiente paso en qiime (feature-table), y este paso a su vez también es necesario para la generación del árbol filogenético. Además, *table.qza* también es necesario para *phyloseq*. :::success **Todo se ejecutó correctamente** :smiley_cat: ::: * **VISUALIZACIÓN DE LOS RESULTADOS TRAS EJECUTAR DADA2** Tras dejarlo corriendo todo el fin de semana, DADA2 se ejecutó correctamente :slightly_smiling_face:, obteniendo finalmente el archivo *table.qza* necesario para los pasos sucesivos. Hecho esto, vamos a explorar los datos obtenidos, usando los comandos siguientes: ``` qiime feature-table tabulate-seqs \ --i-data rep-seqs.qza \ --o-visualization rep-seqs.qzv qiime metadata tabulate \ --m-input-file denoising-stats.qza \ --o-visualization denoising-stats.qzv qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file metadata.tsv ``` :::success **Todo se ejecutó correctamente** :smiley_cat: ::: Mientras tanto, al ejecutar el comando de la ventana gris anterior, obtenemos lo siguiente (pasando los ficheros .qzv al PC local como se ha mostrado anteriormente en este cuaderno y subiéndolos [aquí](https://view.qiime2.org/)): <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero rep-seqs.qzv</i></b></p> </font> ![](https://i.imgur.com/40kjoMo.png) Este fichero nos muestra una tabla con los IDs de las características o *features* y sus secuencias asociadas ([fuente](https://docs.qiime2.org/2022.2/tutorials/moving-pictures/)). <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero denoising-stats.qzv</i></b></p> </font> ![](https://i.imgur.com/qaD9oD1.png) Este fichero nos enseña, para cada muestra, la variación en el número de secuencias tras los sucesivos pasos de filtrado llevados a cabo anteriormente ([fuente](https://forum.qiime2.org/t/interpreting-denoising-stats/6280/2)). <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero table.qzv</i></b></p> </font> ![](https://i.imgur.com/OriBUYU.png) Este fichero nos da información acerca del número de muestras, el número de taxones existentes, el número de taxones identificados por muestra, cómo varían las respuestas en función de la variable seleccionada de los metadatos, etc. * **CONSTRUCCIÓN DEL ÁRBOL FILOGENÉTICO** Para la construcción del árbol filogenético, vamos usar la siguiente *pipeline*, que recurre al programa *mafft*: ``` qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences filtered_rep-seqs.qza \ --o-alignment aligned_filtered-rep-seqs.qza \ --o-masked-alignment masked-aligned-rep-seqs.qza \ --o-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza ``` :::success **Todo se ejecutó correctamente** :smiley_cat: ::: Los archivos generados en este paso serán necesarios para los comandos de *phyloseq*, concretamente el del árbol enraizado (*rooted-tree.qza*). **Nótese, además, que se ha tomado como *input* el fichero *filtered_rep-seqs.qza*, el cual ya ha sido filtrado respecto a las secuencias de mitocondrias y cloroplastos**. * **ASIGNACIÓN TAXONÓMICA** A continuación, vamos a llevar a cabo la asignación taxonómica con SILVA, lo cual es necesario hacer antes de los análisis de alfa y beta diversidad ya que se requiere para dichos pasos un fichero llamado *taxonomy.qza* que se genera con SILVA. Para ello, me he descargado este clasificador de la página de qiime2 (es la versión más reciente): ``` wget https://data.qiime2.org/2022.2/common/silva-138-99-nb-classifier.qza ``` Tras esto, ejecutamos qiime2 para la asignación taxonómica (con el fin de concoer la composición taxonómica de las muestras) usando este clasificador de tipo Naïve-Bayes (esto lleva su tiempo :clock1: :hourglass_flowing_sand:): ``` qiime feature-classifier classify-sklearn \ --i-classifier silva-138-99-nb-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv scp lmarcos@172.16.202.45:~/victor/PYM/taxonomy.qzv . ``` Tras ejecutar estos comandos, obtenemos un fichero .qzv que puede ser visualizado [aquí](https://view.qiime2.org/), observándose lo siguiente: ![](https://i.imgur.com/EK4DSuF.png) Tal y como podemos apreciar, a cada taxón identificado se le asigna un valor de confianza. La mayoría de los taxones están expresados a nivel de género, si bien hay algunos en los que se especifica también la especie. Los taxones con mayor confianza están asociados a mitocondrias/cloroplastos, los cuales han aparecido aquí al poseer estos orgánulos también rRNA 16S. Tal y como discutimos en el GM, **vamos a eliminar estos taxones** que no son de relevancia aquí: ``` qiime taxa filter-table \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --p-exclude mitochondria,chloroplast \ --o-filtered-table filtered_table.qza qiime taxa filter-seqs \ --i-sequences rep-seqs.qza \ --i-taxonomy taxonomy.qza \ --p-exclude mitochondria,chloroplast \ --o-filtered-sequences filtered_rep-seqs.qza ``` Filtrada la **tabla** original, así como las **secuencias representativas**, podemos pedir a qiime que nos dibuje un barplot donde se representen las abundancias de cada ASV para cada caso (aquí he tenido que usar otro fichero llamado metadata.tsv (los datos separados por tabulaciones)): ``` qiime taxa barplot --i-table filtered_table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file datos_IAS.tsv \ --o-visualization taxa_barplot.qzv ``` Si lo visualizamos, obtenemos esto a nivel general, si bien podemos modificar los parámetros según nuestros intereses: <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero taxa_barplot.qzv</i></b></p> </font> ![](https://i.imgur.com/ZzAbVX9.png) Tras obtener el fichero *taxonomy.qza*, así como *filtered_table.qza* y *filtered_rep-seqs.qza*, podemos realizar el análisis de alfa/beta-diversidad en *phyloseq*. * **ANÁLISIS DE LAS ABUNDANCIAS RELATIVAS** A continuación, vamos a llevar a cabo el análisis de las abundancias relativas, así como de las alfa y beta diversidades por medio del paquete *phyloseq* de R. Para ello, ejecutamos R en la consola de Ubuntu e instalamos el paquete, además de otros paquetes necesarios para importar los datos: ``` R #ejecutar R if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("phyloseq") #instalación de phyloseq if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("microbiome") #instalación de microbiome install.packages("remotes") remotes::install_github("jbisanz/qiime2R") #instalación de qiime2R install.packages("readr") #instalación de readr install.packages("tidyr") #instalación de tidyr install.packages("dplyr") #instalación de dplyr install.packages('knitr', dependencies = TRUE) #instalación de knitr conda install -c r r-nloptr install.packages("car") #instalación de car ``` Tras 15 minutos instalándose aproximadamente todos los paquetes, importamos las librerías: ``` library(phyloseq) library(tibble) library(qiime2R) library(tidyr) library(readr) library(dplyr) library(microbiome) library(knitr) library(car) library(stats) ``` Tras esto, pasamos los datos a la terminal de R para posteriormente poder hacer los análisis. Nótese el empleo de la **tabla ya filtrada** para generar la variable *ASV*: ``` #Tabla de OTUs/ASV ASV<-read_qza("./filtered_table.qza") #Taxonomía taxonomyq <-read_qza("./taxonomy.qza") #Transformar tabla 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) #Árbol filogenético que generamos anteriormente tree<-read_qza("./rooted-tree.qza") #Metadata metadata<-read_csv("./metadata.csv") #Construir Objeto phyloseq 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$ID) physeq<- merge_phyloseq(OTUs,tree,TAXq,sample_metadata) ``` Tras esto, he obtenido un objeto llamado *physeq* con el siguiente aspecto: ![](https://i.imgur.com/qbVexdP.png) Los distintos componentes del objeto *physeq* pueden ser consultados de la siguiente manera: ``` otu_table(physeq) #tabla de OTUs sample_data(physeq) #número de muestras tax_table(physeq) #lista de taxones identificados phy_tree(physeq) #información del árbol filogenético nsamples(physeq) #número de muestras sample_names(physeq) #nombre de las muestras taxa_names(physeq) #nombre de los taxones identificados sample_variables(physeq) #nombre de las variables ``` Creado este objeto, antes de analizar las alfa y beta diversidades, vamos a echar un vistazo a cómo las abundancias cambian en función de las distintas variables (o features) de los metadatos. * **ABUNDANCIAS** * **NORMALIZACIÓN DE LOS *COUNTS*** En primer lugar, antes de determinar las abundancias de algunas variables, es necesario que exista una equitativa profundidad de secuenciación entre las muestras, ya que en algunas se han identificado más *features* que en otras (como se puede ver en *table.qzv*). Para conseguir esto, vamos a ejecutar el siguiente comando en R: ``` pseq.rarified <- rarefy_even_depth(physeq) ``` Sin embargo, si ejecutamos *rarefy_even_depth* de esta forma, esta función de R tomará como profundidad, por defecto, el mínimo de features existente entre todas las muestras (en este caso, según *table.qzv*, la muestra PYM307 tiene 47046 *counts* y es la que menos tiene, luego *rarefy_even_depth* tomará dicho número y hará que todas las muestras tengan dicho valor de *counts*) ([FUENTE](https://rdrr.io/bioc/phyloseq/man/rarefy_even_depth.html)). Para saber realmente qué número de counts deberían tener todas las muestras tras esta normalización con *rarefy_even_depth*, vamos a ejectuar el siguiente comando en qiime: ``` qiime diversity alpha-rarefaction --i-table filtered_table.qza \ --p-max-depth 79420 \ --p-steps 20 \ --i-phylogeny rooted-tree.qza \ --m-metadata-file metadata.tsv \ --o-visualization rarefaction_curves.qzv ``` Esta función de qiime2 toma la tabla de datos (ya filtrada) y, sabiendo el número máximo de counts existente, los metadatos y el árbol enraizado, nos permite conocer, para cada variable de los metadatos, qué valor de profundidad de secuenciación es óptima para tener la mayor cantidad de información posible, algo especialmente útil para los análisis de diversidad. Este valor óptimo se conoce en base al punto en el que las gráficas alcanzan pendiente=0. Dicho esto, vamos a ver los resultados para las variables "Nivel_actividad_fisica_total" y "Alcohol", que he seleccionado al existir un balance de muestras para las distintas opciones de dichas variables: * **Nivel de actividad física total** ![](https://i.imgur.com/xnaDOR0.png) * **Alcohol** ![](https://i.imgur.com/0xldUS4.png) Analizando ambas gráficas, podemos ver que, en torno a los **50000** ***counts***, las gráficas se estabilizan para ambas variables y todas sus posibles opciones, luego vamos a tomar, en principio, este valor para realizar la normalización de la profundidad de secuenciación con *rarefy_even_depth*: ``` pseq.rarified <- rarefy_even_depth(physeq, sample.size=50000) ``` :::warning Sin embargo, ejecutar este comando provoca la **pérdida** de aquellas muestras cuyo número de counts sea menor a 50000, algo no deseable. ![](https://i.imgur.com/tLlbGQF.png) ::: Por tanto, para evitar esa pérdida de información, he optado por ejecutar el comando con las **opciones por defecto**. ``` pseq.rarified <- rarefy_even_depth(physeq) ``` * **TRANSFORMACIÓN DE LAS ABUNDANCIAS EN ABUNDANCIAS RELATIVAS** A continuación, tras normalizar el número de *counts* de todas las muestras (a 47406), antes de analizar cómo varían las abundancias taxonómicas en función de las variables, vamos a determinar las **abundancias relativas** (*counts* de un taxón concreto en una muestra / nº total de *counts* para dicha muestra). Para ello, usamos este comando de R: ``` pseq_rel<-transform_sample_counts(pseq.rarified, function (x) x / sum(x) ) ``` Podemos comprobar que las abundancias son ahora relativas de la siguiente forma: ``` head(abundances(pseq.rarified)) head(abundances(pseq_rel)) ``` ![](https://i.imgur.com/7Q24tl8.png) ![](https://i.imgur.com/l7RdhJa.png) Para la muestra **PYM308**, por ejemplo, para el feature **fe6426a2cf4fa1596e6cfec556790cfb**, la abundancia absoluta era antes 23 y ahora es 0.0004888832, que efectivamente es 23/47046, ya que ahora el nº de *counts* de PYM308, así como el de todas las demás muestras, es 47046 tras la normalización con *rarefy_even_depth*. Por tanto, ya hemos normalizado las muestras y calculado las abundancias relativas, englobadas en el objeto *pseq_rel*. Tras esto, vamos a analizar cómo varían estas abundancias en base a distintas variables de los metadatos. Por ejemplo, vamos a determinar cómo varían las abundancias en función del nivel de actividad física de los 70 individuos tenidos en cuenta para este estudio: ``` png(file="Act_fisica.png",width=2000, height=2000) plot_bar(pseq_rel,"Nivel_actividad_fisica_total","Abundance","Phylum",title="Abundancias relativas según el nivel de actividad física total") dev.off() ``` * [**Abundancias relativas según actividad física.**](https://drive.google.com/file/d/1aAspXYio3DwLnCYHiqkA5rYO8AMuIqK6/view?usp=sharing) Tal y como se puede apreciar, la abundancia taxonómica va decreciendo a medida que también disminuye la actividad física del individuo. Por otro lado, vamos a estudiar cómo varían las abundancias en base al consumo de alcohol: ``` png(file="Alcohol.png",width=2000, height=2000) plot_bar(pseq_rel,"Alcohol","Abundance","Phylum",title="Abundancias relativas según el consumo de alcohol") dev.off() ``` * [**Abundancias relativas según el consumo de alcohol.**](https://drive.google.com/file/d/1c7IkUCHizLbme3OoshepTUd0V8x3jqx_/view?usp=sharing) Además, también podemos representar cómo varían los distintos taxones identificados en función de estas variables y los distintos valores de las mismas. ``` plot_bar(pseq_rel,"Phylum","Abundance","Phylum",title="Título",facet_grid="Nivel_actividad_fisica_total~.") plot_bar(pseq_rel,"Phylum","Abundance","Phylum",title="Título",facet_grid="Alcohol~.") plot_bar(pseq_rel,"Phylum","Abundance","Phylum",title="Título",facet_grid="IMC_clasificacion~.") ``` Para el nivel de actividad física o el consumo de alcohol, obtendríamos estas gráficas (he puesto unos enlaces para poder ampliar las imágenes, ya que aquí no se puede): * [**Variación de las abundancias taxonómicas para cada taxón en base a la actividad física.**](https://drive.google.com/file/d/1ZI-QCe0JBU7fbh4gN-rGQeH_W1N6ID4o/view?usp=sharing) * [**Variación de las abundancias taxonómicas para cada taxón en base al consumo de alcohol.**](https://drive.google.com/file/d/19u5gwhZJ9FPaRm55UgcgXe84-z_zLD-P/view?usp=sharing) * [**Variación de las abundancias taxonómicas para cada taxón en base al IMC.**](https://drive.google.com/file/d/1ccoyZGpreXUQrjegI9yn9GohX4tesS1n/view?usp=sharing) :::danger **Traté de ejecutar esto para ignorar, por ejemplo, Firmicutes (porque es el más abundante) en las gráficas anteriores pero las gráficas me salen iguales :crying_cat_face:** ``` filtered<-subset_taxa(pseq_rel,Phylum!="Firmicutes") ``` ::: :::warning **Hay que ver qué hacer con los NA del IMC, por ejemplo, y ver cómo quitar los Phylum más abundantes.** **En el IMC, hay 45 personas con normopeso y 23 con sobrepeso (realmente es sobrepeso+obesidad).** ::: ___ **ANÁLISIS DE ALFA Y BETA DIVERSIDAD A NIVEL GENERAL** Tras esto, vamos a determinar la **alfa-diversidad**. * **ALFA-DIVERSIDAD** ``` rich <- alpha(physeq, index = "all") kable(head(rich)) ``` Este comando genera una tabla con todos los indicadores globales para cada muestra (medidas de riqueza, diversidad, dominancia, rareza, etc., con los parámetros por defecto). Adicionalmente, podemos dibujar gráficas de alfa-diversidad usando distintas medidas de la misma, por medio de la función *plot_richness*: ``` plot_richness(physeq) ``` * [Alfa-diversidad usando distintas medidas](https://drive.google.com/file/d/1kLJG8mFE2tCoqtqWHRjf4zczWjnKA9YG/view?usp=sharing) También podríamos, por ejemplo, evaluar la alfa-diversidad de los individuos en base a una de las variables, como el **nivel de actividad física**, y además analizando solamente ciertas métricas de alfa-diversidad, tal y como se muestra a continuación para "Chao1", "Shannon" y "Simpson": ``` png(file="act_fisica_alpha.png",width=500, height=500) plot_richness(physeq,x="Nivel_actividad_fisica_total",measures=c("Chao1","Shannon","Simpson")) dev.off() ``` ![](https://i.imgur.com/uFQyb1t.png) :::warning **Esperaba que hubiese mayores diferencias entre alto-moderado-bajo, pero no veo diferencias significativas a simple vista.** ::: A continuación, vamos a normalizar el número de *feature counts* para todas las muestras, por medio de *rarefy_even_depth* de *phyloseq*: ``` pseq.rarified <- rarefy_even_depth(physeq) ``` Y, posteriormente, representar la gráfica anterior, por ejemplo, con ese nuevo objeto llamado *pseq.rarified*: ``` png(file="act_fisica_alpha.png",width=500, height=500) plot_richness(pseq.rarified,x="Nivel_actividad_fisica_total",measures=c("Chao1","Shannon","Simpson")) dev.off() ``` ![](https://i.imgur.com/1xEjddM.png) Como vemos, tras la normalización, las gráficas son prácticamente iguales. También podemos representar esta gráfica de la siguiente forma: ``` png(file="act_fisica_alpha_coloured.png",width=500, height=500) plot_richness(pseq.rarified,x="Nivel_actividad_fisica_total",color="Nivel_actividad_fisica_total",measures=c("Chao1","Shannon","Simpson"))+ geom_boxplot(aes(fill = Nivel_actividad_fisica_total), alpha=.7) + theme_bw() dev.off() ``` ![](https://i.imgur.com/LZqJSNi.png) Además, los puntos representados en las gráficas del enlace anterior ([Alfa-diversidad usando distintas medidas](https://drive.google.com/file/d/1kLJG8mFE2tCoqtqWHRjf4zczWjnKA9YG/view?usp=sharing)) (sobre el objeto ya normalizado) pueden ser consultados en una tabla, por medio de: ``` alpha <- estimate_richness(pseq.rarified, measures=c("Chao1", "Shannon", "Simpson")) ``` ![](https://i.imgur.com/tonLGZi.png) Podemos fusionar esta tabla con la tabla de metadatos de la siguiente forma: ``` alpha<- alpha %>% as.data.frame() %>% rownames_to_column("SampleID") %>% left_join(metadata,by=c("SampleID"="ID")) head(alpha) ``` La tabla obtenida tiene el siguiente aspecto (parcialmente): ![](https://i.imgur.com/zXxn3wk.png) Tras generar esta tabla, vamos a determinar si existen diferencias entre las varianzas usando la medida Chao1, por medio del **test de Levene**, y también para la variable que hemos estado analizando, el nivel de actividad física total: ``` leveneTest(y = alpha$Chao1, group = alpha$Nivel_actividad_fisica_total, center = "median") ``` ![](https://i.imgur.com/ICmZAyg.png) Según este test, **no** existen diferencias significativas entre las varianzas de los grupos existentes para el nivel de actividad física total (Alto, Moderado y Bajo). Asimismo, también para esta variable de los metadatos, podemos realizar un **test de Krustal-Wallis**, para determinar si existen diferencias significativas entre los valores de alfa-diversidad, usando Chao1, de los 3 grupos mencionados anteriormente: ``` Chao1 <- kruskal.test(Chao1 ~ Nivel_actividad_fisica_total, data = alpha) ``` ![](https://i.imgur.com/5Fsewq3.png) Tampoco parecen existir diferencias entre las alfa-diversidades de los tres grupos. --- * **BETA-DIVERSIDAD** Transformamos los datos: ``` pseq.compositional <- transform(physeq, "compositional") ``` Y determinamos la beta-diversidad, coloreando también en base a la variable "Nivel_actividad_fisica_total": ``` ordinated_taxa_unifrac = ordinate(pseq.compositional, method="MDS", distance="unifrac",weighted=TRUE) plot_ordination(pseq.compositional, ordinated_taxa_unifrac, color="Nivel_actividad_fisica_total", title = "Unifrac MDS Analysis") + theme_bw() ``` ![](https://i.imgur.com/92CqMW2.png) ___ <font color = 'steelblue'> <p style="text-align:center;"><b>ANÁLISIS DE ALFA Y BETA DIVERSIDAD CON LA TABLA DE GRUPOS SEGÚN CONSUMO DE ULTRAPROCESADOS</b></p> </font> <font color = 'darkorange'> <p style="text-align:center;"><b>Información previa</b></p> </font> :::warning :memo: Nuevos datos a incluir en metadata - Descargar fichero :arrow_down: [aquí ](https://drive.google.com/file/d/1LSMKvswD6gJJres4-8DbIPl_2sNtmODR/view?usp=sharing) **Interpretación Índice de Alimentación Saludable (IAS):** * Dieta excelente: >80 * Muy buena: 71-80 * Buena: 61-70 * Aceptable: 51-60 * Inadecuada <=50 **Cuestionario de Alimentos Altamente Procesados (sQ-HPF)** Equivalencia entre la puntuación Cuestionario de Alimentos Altamente Procesados (sQ-HPF) y el porcentaje estimado de consumo de alimentos altamente procesados sobre la ingesta total en gramos al día. | Puntuación | Porcentaje | | -------- | -------- | | 1 | 11.3 | | 2 | 15 | | 3 | 18.7 | | 4 | 22.4 | | 5 | 26.1 | | 6 | 29.8 | | 7 | 33.5 | | 8 | 37.2 | | 9 | 40.9 | | 10 | 44.6 | | 11 | 48.3 | | 12 | 52 | | 13 | 55.7 | | 14 | 59.4 | ::: :::info :memo: Nota por Laura Hay que limpiar la tabla de taxonomía, quitar los reads que no se han assignado a nada, secuencias asignadas a Eukaryota u otros taxa que sabemos que no deberían estar y que son errores de asignación. Luego aglomeramos al nivel de género, **que es con lo que finalmente vamos a trabajar.** Y transformamos a Abundancia relativas con el paquete `microbiome`. Por último podemos transformar el objeto phyloseq en un dataframe y hacer gráficas chulas con `ggplot` ```r=1 physeq_fil<- subset_taxa(physeq, Kingdom!="Unassigned" & Kingdom!="Eukaryota") phy_gen <- tax_glom(physeq_fil,taxrank="Genus") phy_gen_comp <- microbiome::transform(phy_gen, "compositional") tb_ASV <- psmelt(phy_gen_comp) ``` [Aquí ](https://microbiome.github.io/tutorials/cleaning_taxonomy_table.html) hay algunos tips que pueden servir para limpiar la tabla de taxa :wink: ::: Tras haber generado la tabla separando a los individuos del estudio en dos grupos, y añadiendo información relativa al consumo de alcohol, tabaco, índice de alimentación saludable, etc., vamos a generar un nuevo objeto de *phyloseq* para realizar los análisis de diversidad. ```r=1 # Cargamos las librerías necesarias library(phyloseq) library(tibble) library(qiime2R) library(tidyr) library(readr) library(dplyr) library(microbiome) library(knitr) library(car) library(stats) ``` Tras cargar las librerías, generamos el objeto de *phyloseq*. Para ello, se recurrirá a la tabla de ASVs ya filtrada (es decir, tras haber eliminado lecturas relacionadas con mitocondrias y cloroplastos), la información taxonómica, el árbol enraizado y los datos de los grupos de consumo de ultraprocesados. ```r=1 #Generamos el objeto phyloseq #Tabla de OTUs/ASV ASV<-read_qza("./filtered_table.qza") #Taxonomía taxonomyq <-read_qza("./taxonomy.qza") #Transformar tabla 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) #Árbol filogenético que generamos anteriormente tree<-read_qza("./rooted-tree.qza") #Metadata metadata<-read_csv("./datos_ias.csv") #Construir Objeto phyloseq 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()) sample_metadata <- sample_data(metadata %>% as.data.frame()) sample_names(sample_metadata) <- paste(metadata$SampleID) physeq_processed_foods<- merge_phyloseq(OTUs,tree,TAXq,sample_metadata) ``` **Objeto *physeq_processed_foods***: ![](https://i.imgur.com/FfMIKwL.png) Ahora, vamos a limpiar la tabla de taxonomía, quitando los *reads* que no se han asignado a nada y las secuencias asignadas a *Eukaryota*. Luego aglomeramos al nivel de **género** y transformamos a abundancias relativas con el paquete *microbiome*. ```r=1 #Eliminamos reads no asignados a nada o a Eukaryota. physeq_fil<- subset_taxa(physeq_processed_foods, Kingdom!="Unassigned" & Kingdom!="Eukaryota") #Aglomeramos todo a nivel de género, eliminando además los NAs. phy_gen <- tax_glom(physeq_fil,taxrank="Genus",NArm=TRUE) #Abundancias relativas. phy_gen_comp <- microbiome::transform(phy_gen, "compositional") ``` Tras la eliminación de los taxones indicados anteriormente, pocos taxones son desechados, mientras que dicho número decrece mucho más tras la aglomeración a nivel de género: **Objeto *phy_gen***: ![](https://i.imgur.com/6WWlcFc.png) A continuación, vamos a analizar las alfa-diversidades para los dos grupos de individuos en base al consumo de ultraprocesados. Recordemos que: * **Grupo 1**: individuos cuyo porcentaje estimado de consumo de alimentos altamente procesados sobre la ingesta total es **≤15%**. * **Grupo 2**: individuos cuyo porcentaje estimado de consumo de alimentos altamente procesados sobre la ingesta total es **>15%**. ```r=1 library(dplyr) library(tidyr) library(ggplot2) library(ggpubr) library(janitor) library(readxl) library(phyloseq) library(ape) library(mia) library(vegan) library(rstatix) library(scater) library(ggsignif) p_alpha <- makeTreeSummarizedExperimentFromPhyloseq(phy_gen) tse.list <- list("p_alpha" = p_alpha) tse.list <- lapply(names(tse.list), function(x){ tse.list[[x]] <- tse.list[[x]] %>% estimateRichness(abund_values = "counts", index = "observed", name = "observed") %>% estimateRichness(abund_values = "counts", index = "chao1", name = "chao1") %>% estimateDiversity(abund_values = "counts", index = "shannon", name = "shannon") %>% estimateDiversity(abund_values = "counts", index = "gini_simpson", name = "simpson") %>% estimateEvenness(abund_values = "counts", index = "pielou", name = "pielou") %>% estimateDominance(abund_values = "counts", index = "relative", name = "relative") %>% estimateDiversity(abund_values = "counts", index = "log_modulo_skewness", name = "log_modulo_skewness") }) names(tse.list) <- c("PYM") df_alpha <- as.data.frame(colData(tse.list$PYM)) %>% select(chao1, shannon, simpson, pielou, relative, log_modulo_skewness, Grupo_HPF) counter <<- 0 graphs_alpha <-lapply(df_alpha[ ,c("chao1", "shannon", "simpson", "pielou", "relative", "log_modulo_skewness")], function(a) {counter <<- counter + 1 ggplot(df_alpha, aes(x = Grupo_HPF, y = a)) + geom_boxplot(aes(fill = Grupo_HPF), alpha=.5, outlier.shape = NA) + geom_signif(comparisons = list(c("Group 1", "Group 2")), test = "wilcox.test", map_signif_level = TRUE, textsize = 3, fontface = "bold") + geom_jitter(width = 0.2, aes(colour = Grupo_HPF), size = 1.5) + ylab(gsub('_', ' ', colnames(df_alpha)[counter])) + xlab(NULL) + theme_bw() + theme(axis.title.y = element_text(size=12, face="bold.italic", colour = "black"), axis.text.x = element_text(size=9, face="bold", colour = "black"), axis.text.y = element_text(size=9, face="bold", colour = "black"), legend.title=element_blank(), legend.text = element_text(face = "bold"))}) ggarrange(plotlist = graphs_alpha, common.legend = TRUE, legend = "right", align = "hv") ``` :::warning **Las alfa-diversidades se hacen con las abundancias absolutas** ::: <font color = 'black'> <p style="text-align:center;"><b>ALFA-DIVERSIDADES</b></p> </font> ![](https://i.imgur.com/nlCzQC1.png) Como podemos apreciar, parece existir algo más de homogeneidad entre las alfa-diversidades de las muestras del grupo 1 en comparación con el grupo 2, si bien no parece haber diferencias significativas entre ambos grupos. :::warning **CORRECCIÓN POR *MULTIPLE TESTING*:** ```r=1 wil.test_alpha <- bind_rows(pairwise_wilcox_test(df_alpha, chao1 ~ Grupo_HPF), pairwise_wilcox_test(df_alpha, shannon ~ Grupo_HPF), pairwise_wilcox_test(df_alpha, simpson ~ Grupo_HPF), pairwise_wilcox_test(df_alpha, pielou ~ Grupo_HPF), pairwise_wilcox_test(df_alpha, relative ~ Grupo_HPF), pairwise_wilcox_test(df_alpha, log_modulo_skewness ~ Grupo_HPF)) %>% adjust_pvalue(method = "BH") %>% add_significance() ``` ``` wil.test_alpha # A tibble: 6 x 9 .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr> 1 chao1 Group 1 Group 2 27 33 480. 0.618 0.742 ns 2 shannon Group 1 Group 2 27 33 503 0.4 0.626 ns 3 simpson Group 1 Group 2 27 33 501 0.417 0.626 ns 4 pielou Group 1 Group 2 27 33 512 0.329 0.626 ns 5 relative Group 1 Group 2 27 33 365 0.236 0.626 ns 6 log_modulo_skewness Group 1 Group 2 27 33 435 0.883 0.883 ns ``` **Nada significativo, luego la gráfica de antes está bien :smile_cat:** ::: Esto lo vamos a comprobar también mediante el análisis de la **beta-diversidad**: ```r=1 library(dplyr) library(tidyr) library(ggplot2) library(ggpubr) library(RColorBrewer) library(janitor) library(readxl) library(ape) library(mia) library(vegan) library(rstatix) library(scater) library(phyloseq) library(tibble) library(qiime2R) library(readr) library(microbiome) library(knitr) library(car) library(stats) library(cowplot) distances <- function(study_pseq) { for( i in dist_methods ){ # Calculate distance matrix iDist <- phyloseq::distance(study_pseq, method = i) dlist[[i]] <- iDist # Calculate ordination iPCoA <- ordinate(study_pseq, "PCoA", distance = iDist) pcoa_list[[i]] <- iPCoA } to_return <- list('dlist' = dlist, 'pcoa_list' = pcoa_list) return(to_return) } # Performs PERMANOVA: adonis_calculator <- function(dlist, study_pseq) { results.group <- lapply(names(dlist), function(x) { z <- adonis2(dlist[[x]] ~ Grupo_HPF, data = data.frame(sample_data(study_pseq))) return(as.data.frame(z)) }) names(results.group) <- names(dlist) return(results.group) } # Makes plots: plotter_beta <- function(dist_methods, study_pseq, pcoa_list){ for( i in dist_methods ){ p <- NULL p <- plot_ordination(study_pseq, pcoa_list[[i]], color = "Grupo_HPF") + geom_point(size = 2) + stat_ellipse(aes(group = Grupo_HPF), linetype = 2) + theme_bw() + theme(plot.title = element_text(face = 'bold', size = 16), axis.title = element_text(size = 14), legend.title = element_blank(), axis.text = element_text(size = 12), legend.text = element_text(size = 12, face = "bold")) p <- p + ggtitle(i) p <- p + scale_colour_brewer(type="qual", palette="Set2") plist[[i]] <- p } return(plist) } dist_methods <- unlist(distanceMethodList)[c(1, 2, 8, 10)] dist_methods # to do this individually for each study: dlist <- vector("list", length(dist_methods)) # distance matrix names(dlist) <- dist_methods pcoa_list <- dlist # PCoA plist <- dlist # plots phy_gen_comp <- microbiome::transform(phy_gen, "compositional") phy_gen_comp.beta <- distances(phy_gen_comp) phy_gen_comp.beta.adonis <- adonis_calculator(phy_gen_comp.beta$dlist, phy_gen_comp) phy_gen_comp.beta.plot <- plotter_beta(dist_methods, phy_gen_comp, phy_gen_comp.beta$pcoa_list) plist <- plotter_beta(dist_methods, phy_gen_comp, phy_gen_comp.beta$pcoa_list) ggarrange(plotlist = plist, common.legend = TRUE, legend = "right") ggarrange(plist[[1]] + ggtitle('a'), plist[[2]] + ggtitle('b'), plist[[3]] + ggtitle('c'), plist[[4]] + ggtitle('d'), common.legend = TRUE, legend = "bottom", align = 'hv') ``` :::warning **La beta-diversidad se ha determinado en base al objeto *phyloseq* con los taxones no asignados ya eliminados, tras aglomerar a nivel de género y tras transformar las abundancias en relativas** ::: <font color = 'black'> <p style="text-align:center;"><b>BETA-DIVERSIDAD</b></p> </font> ![](https://i.imgur.com/iNlj3BH.png) Como podemos apreciar, parece ser que no existen diferencias significativas en la composición del microbioma entre ambos grupos, al menos para la mayoría de las muestras. Vamos a comprobarlo con el **PERMANOVA**: ``` > phy_gen_comp.beta.adonis $unifrac Df SumOfSqs R2 F Pr(>F) Grupo_HPF 1 0.07260026 0.01891883 1.118452 0.317 Residual 58 3.76486080 0.98108117 NA NA Total 59 3.83746106 1.00000000 NA NA $wunifrac Df SumOfSqs R2 F Pr(>F) Grupo_HPF 1 0.03902501 0.01554208 0.9156719 0.419 Residual 58 2.47190124 0.98445792 NA NA Total 59 2.51092624 1.00000000 NA NA $bray Df SumOfSqs R2 F Pr(>F) Grupo_HPF 1 0.08551114 0.01473742 0.8675557 0.516 Residual 58 5.71680429 0.98526258 NA NA Total 59 5.80231543 1.00000000 NA NA $jaccard Df SumOfSqs R2 F Pr(>F) Grupo_HPF 1 0.1688044 0.01589138 0.9365838 0.445 Residual 58 10.4535840 0.98410862 NA NA Total 59 10.6223884 1.00000000 NA NA ``` :::info **Nada significativo, como podía esperarse (mirando los p.valores y R2).** ::: --- Además, podemos determinar qué géneros son más abundantes en los dos grupos y ver si hay diferencias en las abundancias entre dichos grupos, por medio de: ```r=1 #Porcentaje Genus tb_gen <- psmelt(phy_gen_comp) tbgen <- tb_gen %>% group_by(Genus) %>% summarise(M=median(Abundance)) %>% ungroup() %>% unique() %>% top_n(M,n=7) tbgen <- tbgen %>% inner_join(tb_gen) Plotf <- ggplot(tbgen, aes(Grupo_HPF, Abundance ,fill=Genus)) + geom_col(position="fill") + scale_fill_manual(values=c("red","green","blue","orange","yellow","pink","brown")) Plotf <- Plotf + theme_minimal() + xlab("Grupo_HPF") + 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 ``` <font color = 'black'> <p style="text-align:center;"><b>RESULTADO</b></p> </font> ![](https://i.imgur.com/cKw4WQ2.png) En esta gráfica, se están representando los 7 géneros más abundantes (en base a la mediana de sus abundancias relativas para todas las muestras) para cada grupo. Se puede apreciar que, en general, esta mediana de abundancias relativas es similar para casi todos los géneros entre los dos grupos, aunque parece haber algo más de abundancia de *Prevotella* para el grupo 2. Los géneros están apilados por orden alfabético y el tamaño de los bloques en cada barra es relativo a la desviación de las abundancias, en cada grupo, respecto a la mediana de abundancia de dicho género en tal grupo. ---