TFM Víctor - HackMD

# TFM Víctor ###### tags: `Estudiantes` `ML` `TFM` ## **Instalación de paquetes y herramientas** ``` ssh estudiantes@172.16.235.157 #meto la password ``` He instalado conda y he creado el siguiente environment: ``` conda create --name metaphlan ``` Tras esto, he instalado MetaPhlAn (versión 4) y Humann (versión 3) en dicho environment de conda: ``` conda activate metaphlan conda install -c bioconda metaphlan=4.0.0 #OJO: hay que tener una versión de Python > 3.7 ``` Para la instalación de Humann 3.0 (en el mismo environment): ``` conda config --add channels defaults conda config --add channels bioconda conda config --add channels conda-forge conda config --add channels biobakery conda install humann -c biobakery ``` Test de que todo está bien: ``` humann_test #OK humann --version #humann v3.5 metaphlan --version #MetaPhlAn version 4.0.0 (22 Aug 2022) ``` Las versiones de Humann, MetaPhlAn y StrainPhlAn son las más recientes. ## **1ª tanda muestras AI4Food (info)** ### 29/09/2022 Para prácticamente todas las muestras, se aprecia una calidad de secuenciación muy favorable (Phred score: 36 aprox), fiabilidad del 99.99% (sección de data quality control). Ejemplo: ![](https://i.imgur.com/XBNB1C9.png) :::warning ¿Por qué se produce ese salto de calidad a mitad de las lecturas? ::: Por otro lado, respecto a la tasa de error de secuenciación, se observa, para una de las muestras tomada a modo de ejemplo, lo siguiente: ![](https://i.imgur.com/Sggg3kk.png) Según indica el propio informe, la tasa de error crece a medida que se avanza en la secuenciación debido al consumo del reactivo de dicho proceso. :::warning Pasa algo similar a lo de antes, la tasa de error crece hasta llegar a la mitad, y entonces baja bruscamente. ::: Algo similar ocurre también a la hora de determinar la distribución de bases A/G/C/T: ![](https://i.imgur.com/qE50HdG.png) :::info Toda esta información indica que se trata de **secuencias paired-end**. ::: Para todas estas muestras, encontramos una calidad de secuenciación muy alta, en todos los casos superior al 99%. También tenemos los adaptadores: ![](https://i.imgur.com/ov7F0MY.png) También se nos informa, para cada muestra, sobre el número de lecturas, el %GC, la tasa de error de secuenciación (en todos los casos 0.03), etc. --- ## **PRUEBA CON DOS MUESTRAS** A la espera de las muestras de verdad, voy a ir practicando los pasos a realizar con otras muestras de prueba (son 4: **D1800_1.fastq, D1800_2.fastq, D1801_1.fastq, D1801_2.fastq**). :::info **Copiado de Laura** ### **Análisis metagenomas WGS** Yo he estado trabajando con [MetaWRAP](https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md). Mi idea era seguir este pipeline completo, pero no tenemos capacidad para hacer la asignación de MAGs :sweat:. Así que decidí hacer asignación taxonómica con metaphlan... y ya puedo luego comparar mis datos con lo que esté en *curatedMetagenomicData* :smirk: El procesamiento para *metaphlan* :custard: que he usado en mi proyecto de celiacos es el siguiente: 1. Ver QC y quitar adaptadores.... (MetaWRAP usa Read_QC y Trimgallore) 2. Filtrar seq. humanas. **Con BMTagger** Descargar: hg38 index. 3. Asignación taxonómica - [Metaphlan3](https://github.com/biobakery/MetaPhlAn/wiki/MetaPhlAn-3.0) :custard: 4. Análisis funcional - Humann3 :woman_in_lotus_position: > **Nota sobre paired-end** MetaPhlAn can also natively handle paired-end metagenomes (but does not use the paired-end information), and, more generally, metagenomes stored in multiple files (but you need to specify the --bowtie2out parameter): ::: En primer lugar, vamos a instalar **metaWRAP**, siguiendo las instrucciones del github, descargando e indexando asimismo el hg38. --- ### **PASO 1: METAWRAP** **metaWRAP sirve para realizar controles de calidad antes y después de eliminar lecturas de baja calidad y los adaptadores**, recurriendo a TrimGalore y BMTagger. **Instalación:** ```r=1 #me ha costado la misma vida pero #después de 4000 intentos lo conseguí conda install -y mamba git clone https://github.com/bxlab/metaWRAP.git cd metaWRAP mkdir BMTAGGER_INDEX cd BMTAGGER_INDEX wget ftp://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/*fa.gz gunzip *fa.gz cat *fa > hg38.fa rm chr*.fa bmtool -d hg38.fa -o hg38.bitmask srprism mkindex -i hg38.fa -o hg38.srprism #esto tarda un rato export PATH=$PATH:/home/estudiantes/victor/metaWRAP/bin mamba create -y -n metawrap python=2.7 conda activate metawrap conda config --add channels defaults conda config --add channels conda-forge conda config --add channels bioconda conda config --add channels ursky mamba install --only-deps -c ursky metawrap-mg ``` Una vez instalado metaWRAP, lo ejecutamos tal y como viene indicado en [el tutorial correspondiente](https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md). Sin embargo, al menos con estas dos muestras, los últimos reads en los 4 archivos están corruptos(la longitud de la secuencia es mayor a la secuencia de calidad). Por tanto, para los 4 archivos, es necesario borrar este último read defectuoso para que metaWRAP no dé error (pongo solo un ejemplo para que se vea lo que hice): ```r=1 tail -10 ./RAW_READS/D1800_1.fastq #veo las últimas 10 líneas del fichero head -n -4 ./RAW_READS/D1800_1.fastq > prueba.fastq #muevo todas las lecturas menos la defectuosa a otro archivo ``` Tras hacer lo mismo con los otros 3 archivos .fastq, ejecuto metaWRAP de la siguiente manera: ```r=1 mkdir READ_QC #dentro de la carpeta metaWRAP for F in RAW_READS/PROC_RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 1 -o READ_QC/$SAMPLE & done ``` Tras esto, se generan en la carpeta READ_QC dos carpetas: D1800 y D1801, en donde podemos obtener los FASTQC de las muestras tanto antes como después de ejecutar TrimGalore y BMTagger. Dentro de cada una de las carpetas D1800 y D1801, tenemos las lecturas limpias: ```r=1 final_pure_reads_1.fastq final_pure_reads_2.fastq ``` Estas lecturas las moveremos a una carpeta llamada CLEAN_READS: ```r=1 #dentro de la carpeta metaWRAP mkdir CLEAN_READS for i in READ_QC/*; do b=${i#*/} mv ${i}/final_pure_reads_1.fastq CLEAN_READS/${b}_1.fastq mv ${i}/final_pure_reads_2.fastq CLEAN_READS/${b}_2.fastq done ``` Tras esto, dentro de la carpeta metaWRAP, tenemos una carpeta llamada CLEAN_READS con las lecturas limpias de todas las muestras. --- ### **PASO 2: METAPHLAN4** Ya tenemos las lecturas limpias, ahora vamos a determinar su **composición taxonómica** con **MetaPhlAn4**. **Instalación**: ```r=1 conda create --name mpa -c conda-forge -c bioconda python=3.7 metaphlan conda install -c conda-forge -c bioconda metaphlan conda activate mpa metaphlan --version #MetaPhlAn version 4.0.2 (22 Sep 2022) ``` Tras la instalación, tenemos que instalar también la base de datos **bowtie2**: ```r=1 cd ~/victor metaphlan --install --bowtie2db ~/victor/bowtiedb ``` Tras 1000 horas esperando, todo parece haber ido bien: ```r=1 Download complete The database is installed ``` Ahora sí, vamos a ejecutar MetaPhlAn4 con nuestras muestras de prueba (**nótese que a MetaPhlAn4 se le dan solo como argumento los _1.fastq**): ```r=1 cd ~/victor/metaWRAP/CLEAN_READS touch metaphlan_script.pl #creo un script chmod a+x metaphlan_script.pl #para poder modificarlo vi metaphlan_script.pl #lo abro para modificarlo #copio lo siguiente en dicho fichero: #!/usr/bin/perl #dar de argumento solo los _1 foreach $ar (@ARGV) { # sacamos cada archivo: @field = split (/_/, $ar); $name = $ar; # name = nombre completo $name =~ s/_1.fastq//; # name = nombre antes de _ # obtenemos el _2 correspondiente: $r1 = $ar; # nombre completo de lectura con _1 $r2 = $ar; # lo mismo que r1 por ahora $r2 =~ s/1.fastq/2.fastq/; # cambia _1 por _2 # pasamos ambos a metaphlan: system("metaphlan $r1,$r2 --bowtie2db ~/victor/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt"); } ./metaphlan_script.pl *_1.fastq > metaphlanOUT_26_10_2022_12_30.txt 2>&1 #esto lo hago por si sale algún error. ``` Tras tardar como 20 minutos por muestra, en la carpeta CLEAN_READS, tengo dos ficheros más: **profiled_D1800.txt** y **profiled_D1801.txt**. Los abro con: ```r=1 less profiled_D1800.txt less profiled_D1801.txt ``` Me sale info taxonómica, con abundancias relativas y tal, así que supongo que se ha hecho bien. --- ### **PASO 3: ALFA Y BETA DIVERSIDADES** Para esto, hay que transformar la info obtenida en un objeto de **phyloseq** usando R. --- ## **DATOS AI4FOOD** Ya han llegado las lecturas de las primeras muestras de AI4FOOD, de manera que voy a repetir los pasos previos con todas estas muestras: --- Antes de nada, cada vez que abras la consola tienes que poner: ```r=1 source ~/.bashrc #esto permite que puedas activar los envs ``` --- ### **PASO 1: COPIA CORRECTA DE LOS ARCHIVOS** Antes de nada, tenemos que ver que la copia de los archivos se ha realizado correctamente. TASK: Comprobar que los ficheros se han copiado correctamente Dentro de la ruta victor/AI4FOOD/RAW_READS, tengo todos los ficheros fq.gz, que he copiado desde el disco duro que nos enviaron. Vamos a ejecutar los siguientes comandos, usando md5sum, para saber que la copia se ha hecho bien. Dentro de esta carpeta RAW_READS (que es donde están todos los fq.gz), ejecuto los siguientes comandos: ```r=1 for file in $(find -type f); do echo $file; md5sum $file > $file.md5; done ``` ```r=1 for file in $(find -type f | grep md5); do cat $file >> md5_after_unzipping.txt ; done ``` Compruebo que todo se ha realizado bien: ```r=1 md5sum --check MD5.txt > md5sum_check_log.txt 2>&1 ``` Para todos los casos sale: ``` La suma coincide ``` Todo guay, luego pasamos al paso 2: --- ### **PASO 2: UNIR MUESTRAS SECUENCIADAS DOS VECES** TASK: Unir aquellos ficheros que vengan de la misma muestra y sentido de secuenciación En primer lugar, vamos a descomprimir los ficheros: ```r=1 gunzip *.gz #tarda 1000 horas ``` Tras esto, obtenemos varios ficheros .fq, detectando varias muestras secuenciadas dos veces, las cuales son las siguientes: | Muestra | Detalles | | -------- | -------- | | **V1_614** | Secuenciada 2 veces (L1 y L3, también cambia el flowcell)|| -------- | -------- | | **V1_650** | Secuenciada 2 veces (L2 y L3)|| -------- | -------- | | **V1_668** | Secuenciada 2 veces (L1 y L3, también cambia el flowcell)|| -------- | -------- | | **V1_752** | Secuenciada 2 veces (L1 y L3, también cambia el flowcell)|| -------- | -------- | | **V1_760** | Secuenciada 2 veces (L1 y L3, también cambia el flowcell)|| -------- | -------- | | **V1_783** | Secuenciada 2 veces (L3 y L4, también cambia el flowcell)| | **V1_817** | Secuenciada 2 veces (L1 y L4, también cambia el flowcell)| | **V1_1120** | Secuenciada 2 veces (L1 y L4, también cambia el flowcell)| | **V1_1251** | Secuenciada 2 veces (L2 y L3, también cambia el flowcell)| | **V1_1310** | Secuenciada 2 veces (el lane es el mismo (L1) pero cambia el flowcell (en uno es H3NKVDSX5_L1 y en otro es H5J5KDSX5_L1))| | **V1_1797** | Secuenciada 2 veces (el lane es el mismo (L1) pero cambia el flowcell (en uno es H3NKVDSX5_L1 y en otro es H5J5KDSX5_L1))| | **V1_1976** | Secuenciada 2 veces (el lane es el mismo (L1) pero cambia el flowcell (en uno es H3NKVDSX5_L1 y en otro es H5J5KDSX5_L1))| | **V1_1988** | Secuenciada 2 veces (el lane es el mismo (L1) pero cambia el flowcell (en uno es H3NKVDSX5_L1 y en otro es H5J5KDSX5_L1))| | **V1_1993** | Secuenciada 2 veces (el lane es el mismo (L1) pero cambia el flowcell (en uno es H3NKVDSX5_L1 y en otro es H5J5KDSX5_L1))| | **V1_2131** | Secuenciada 2 veces (L1 y L3, también cambia el flowcell)| | **V1_2157** | Secuenciada 2 veces (L2 y L3, también cambia el flowcell)| | **V1_2167** | Secuenciada 2 veces (el lane es el mismo (L3) pero cambia el flowcell (en uno es H3NKVDSX5_L3 y en otro es H5J5KDSX5_L3))| | **V1_2168** | Secuenciada 2 veces (L3 y L4, también cambia el flowcell)| | **V1_2170** | Secuenciada 2 veces (L2 y L3, también cambia el flowcell)| | **V1_2172** | Secuenciada 2 veces (el lane es el mismo (L3) pero cambia el flowcell (en uno es H3NKVDSX5_L3 y en otro es H5J5KDSX5_L3))| |**V1_2141_a** | Secuenciada 2 veces (el lane es el mismo (L3) pero cambia el flowcell (en uno es H3NKVDSX5_L3 y en otro es H5J5KDSX5_L3))| Una vez identificadas estas muestras secuenciadas más de una vez, tenemos que unir los archivos correspondientes, para así tener, para todas las muestras, dos archivos .fq, y hecho esto ya podremos acortar los nombres y ejecutar metaWRAP. ```r=1 #en la carpeta AI4FOOD mkdir Lecturas_repetidas #tras esto, muevo todos los ficheros de muestras #secuenciadas 2 veces a esta carpeta, usando FileZilla. #Unimos los archivos: cd Lecturas_repetidas cat V1_1120_EKDN220029060-1A_H3NKVDSX5_L1_1.fq V1_1120_EKDN220029060-1A_H5J5KDSX5_L4_1.fq > V1_1120_1.fq cat V1_1120_EKDN220029060-1A_H3NKVDSX5_L1_2.fq V1_1120_EKDN220029060-1A_H5J5KDSX5_L4_2.fq > V1_1120_2.fq cat V1_1251_EKDN220029064-1A_H3YFWDSX5_L2_1.fq V1_1251_EKDN220029064-1A_H5JLKDSX5_L3_1.fq > V1_1251_1.fq cat V1_1251_EKDN220029064-1A_H3YFWDSX5_L2_2.fq V1_1251_EKDN220029064-1A_H5JLKDSX5_L3_2.fq > V1_1251_2.fq cat V1_1310_EKDN220029070-1A_H3NKVDSX5_L1_1.fq V1_1310_EKDN220029070-1A_H5J5KDSX5_L1_1.fq > V1_1310_1.fq cat V1_1310_EKDN220029070-1A_H3NKVDSX5_L1_2.fq V1_1310_EKDN220029070-1A_H5J5KDSX5_L1_2.fq > V1_1310_2.fq cat V1_1797_EKDN220029078-1A_H3NGLDSX5_L1_1.fq V1_1797_EKDN220029078-1A_H5J5KDSX5_L1_1.fq > V1_1797_1.fq cat V1_1797_EKDN220029078-1A_H3NGLDSX5_L1_2.fq V1_1797_EKDN220029078-1A_H5J5KDSX5_L1_2.fq > V1_1797_2.fq cat V1_1976_EKDN220029080-1A_H3NKVDSX5_L1_1.fq V1_1976_EKDN220029080-1A_H5J5KDSX5_L1_1.fq > V1_1976_1.fq cat V1_1976_EKDN220029080-1A_H3NKVDSX5_L1_2.fq V1_1976_EKDN220029080-1A_H5J5KDSX5_L1_2.fq > V1_1976_2.fq cat V1_1988_EKDN220029082-1A_H3NKVDSX5_L1_1.fq V1_1988_EKDN220029082-1A_H5J5KDSX5_L1_1.fq > V1_1988_1.fq cat V1_1988_EKDN220029082-1A_H3NKVDSX5_L1_2.fq V1_1988_EKDN220029082-1A_H5J5KDSX5_L1_2.fq > V1_1988_2.fq cat V1_1993_EKDN220029084-1A_H3NKVDSX5_L1_1.fq V1_1993_EKDN220029084-1A_H5J5KDSX5_L1_1.fq > V1_1993_1.fq cat V1_1993_EKDN220029084-1A_H3NKVDSX5_L1_2.fq V1_1993_EKDN220029084-1A_H5J5KDSX5_L1_2.fq > V1_1993_2.fq cat V1_2131_EKDN220029092-1A_H3NKVDSX5_L1_1.fq V1_2131_EKDN220029092-1A_H5J5KDSX5_L3_1.fq > V1_2131_1.fq cat V1_2131_EKDN220029092-1A_H3NKVDSX5_L1_2.fq V1_2131_EKDN220029092-1A_H5J5KDSX5_L3_2.fq > V1_2131_2.fq cat V1_2141_a_EKDN220029130-1A_H35H7DSX5_L3_1.fq V1_2141_a_EKDN220029130-1A_H5CVJDSX5_L3_1.fq > V1_2141a_1.fq cat V1_2141_a_EKDN220029130-1A_H35H7DSX5_L3_2.fq V1_2141_a_EKDN220029130-1A_H5CVJDSX5_L3_2.fq > V1_2141a_2.fq cat V1_2157_EKDN220029114-1A_H3YFWDSX5_L2_1.fq V1_2157_EKDN220029114-1A_H5JLKDSX5_L3_1.fq > V1_2157_1.fq cat V1_2157_EKDN220029114-1A_H3YFWDSX5_L2_2.fq V1_2157_EKDN220029114-1A_H5JLKDSX5_L3_2.fq > V1_2157_2.fq cat V1_2167_EKDN220029124-1A_H35H7DSX5_L3_1.fq V1_2167_EKDN220029124-1A_H5CVJDSX5_L3_1.fq > V1_2167_1.fq cat V1_2167_EKDN220029124-1A_H35H7DSX5_L3_2.fq V1_2167_EKDN220029124-1A_H5CVJDSX5_L3_2.fq > V1_2167_2.fq cat V1_2168_EKDN220029125-1A_H35H7DSX5_L4_1.fq V1_2168_EKDN220029125-1A_H5CVJDSX5_L3_1.fq > V1_2168_1.fq cat V1_2168_EKDN220029125-1A_H35H7DSX5_L4_2.fq V1_2168_EKDN220029125-1A_H5CVJDSX5_L3_2.fq > V1_2168_2.fq cat V1_2170_EKDN220029127-1A_H35H7DSX5_L3_1.fq V1_2170_EKDN220029127-1A_H3WWHDSX5_L2_1.fq > V1_2170_1.fq cat V1_2170_EKDN220029127-1A_H35H7DSX5_L3_2.fq V1_2170_EKDN220029127-1A_H3WWHDSX5_L2_2.fq > V1_2170_2.fq cat V1_2172_EKDN220029129-1A_H35H7DSX5_L3_1.fq V1_2172_EKDN220029129-1A_H5CVJDSX5_L3_1.fq > V1_2172_1.fq cat V1_2172_EKDN220029129-1A_H35H7DSX5_L3_2.fq V1_2172_EKDN220029129-1A_H5CVJDSX5_L3_2.fq > V1_2172_2.fq cat V1_614_EKDN220029036-1A_H3NKVDSX5_L1_1.fq V1_614_EKDN220029036-1A_H5J5KDSX5_L3_1.fq > V1_614_1.fq cat V1_614_EKDN220029036-1A_H3NKVDSX5_L1_2.fq V1_614_EKDN220029036-1A_H5J5KDSX5_L3_2.fq > V1_614_2.fq cat V1_650_EKDN220029038-1A_H5CW3DSX5_L2_1.fq V1_650_EKDN220029038-1A_H5J5KDSX5_L3_1.fq > V1_650_1.fq cat V1_650_EKDN220029038-1A_H5CW3DSX5_L2_2.fq V1_650_EKDN220029038-1A_H5J5KDSX5_L3_2.fq > V1_650_2.fq cat V1_668_EKDN220029040-1A_H3NKVDSX5_L1_1.fq V1_668_EKDN220029040-1A_H5J5KDSX5_L3_1.fq > V1_668_1.fq cat V1_668_EKDN220029040-1A_H3NKVDSX5_L1_2.fq V1_668_EKDN220029040-1A_H5J5KDSX5_L3_2.fq > V1_668_2.fq cat V1_752_EKDN220029041-1A_H3NKVDSX5_L1_1.fq V1_752_EKDN220029041-1A_H5J5KDSX5_L3_1.fq > V1_752_1.fq cat V1_752_EKDN220029041-1A_H3NKVDSX5_L1_2.fq V1_752_EKDN220029041-1A_H5J5KDSX5_L3_2.fq > V1_752_2.fq cat V1_760_EKDN220029042-1A_H3NKVDSX5_L1_1.fq V1_760_EKDN220029042-1A_H5J5KDSX5_L3_1.fq > V1_760_1.fq cat V1_760_EKDN220029042-1A_H3NKVDSX5_L1_2.fq V1_760_EKDN220029042-1A_H5J5KDSX5_L3_2.fq > V1_760_2.fq cat V1_783_EKDN220029044-1A_H3WWHDSX5_L4_1.fq V1_783_EKDN220029044-1A_H57WKDSX5_L3_1.fq > V1_783_1.fq cat V1_783_EKDN220029044-1A_H3WWHDSX5_L4_2.fq V1_783_EKDN220029044-1A_H57WKDSX5_L3_2.fq > V1_783_2.fq cat V1_817_EKDN220029048-1A_H3NKVDSX5_L1_1.fq V1_817_EKDN220029048-1A_H5J5KDSX5_L4_1.fq > V1_817_1.fq cat V1_817_EKDN220029048-1A_H3NKVDSX5_L1_2.fq V1_817_EKDN220029048-1A_H5J5KDSX5_L4_2.fq > V1_817_2.fq ``` Tras ejecutar esto, obtengo, para cada muestra, dos archivos .fq (uno _1 y otro _2), cuyo tamaño es la suma de los archivos _1 y _2, lo cual indica que el comando ha hecho lo que yo quería. Como vemos, con los comandos anteriores, ya se genera directamente el fichero que yo quiero usar en metaWRAP con el nombre simplificado, lo cual haré también para el resto de muestras que están en la carpeta "RAW_READS". Hecho esto, pues, vamos a pasar al paso 3: --- ### **PASO 3: metaWRAP** TASK: metaWRAP Después de unir los ficheros de muestras que hayan sido secuenciadas más de una vez, vamos a ejecutar metaWRAP. Primero, ejecuto esto para que se pueda ejecutar metaWRAP: ```r=1 export PATH=$PATH:/home/estudiantes/victor/metaWRAP/bin ``` Luego, creo una carpeta READ_QC dentro de AI4FOOD: ```r=1 mkdir READ_QC #dentro de la carpeta AI4FOOD ``` Antes de ejecutar metaWRAP, vamos a acortar los nombres de las muestras para que sea más fácil trabajar con ellas. *Esto **NO** afecta a la integridad de los archivos (he comprobado los md5 y no cambian)*: ```r=1 for f in *.fq; do mv "$f" \ "$(awk -F'_' '{print $1"_"$2"_"$6}' <<< $f)" done ``` Tras acortar los nombres, he comprobado los md5 de algunos archivos, así como su tamaño (usando FileZilla) y todo está en orden, así que el cambio de nombre no ha afectado a la integridad de los archivos. **Al final, todos los ficheros tienen la misma estructura: V1_XXXX_X.fq** (hubo un par de ficheros a los que les tuve que cambiar el nombre manualmente para que el script anterior funcionase). :::info **TEST ANTES DE USAR metaWRAP CON TODAS LAS MUESTRAS** Antes de ejecutar metaWRAP con todas las muestras, voy a hacerlo con un par de muestras (V1_1065 y V1_1100, que están en la carpeta "prueba_metawrap"), para ver si todo va bien: ```r=1 #ACORTO LOS NOMBRES PARA ESAS DOS MUESTRAS cd prueba_metawrap for f in *.fq; do mv "$f" \ "$(awk -F'_' '{print $1"_"$2"_"$6}' <<< $f)" done ``` Ahora tenemos V1_1065_1.fq, V1_1065_2.fq, V1_1100_1.fq y V1_1100_2.fq. Vamos a ejecutar **metaWRAP** para estos dos archivos (desde la carpeta "AI4FOOD"): ```r=1 for F in prueba_metawrap/*_1.fq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 1 -o READ_QC/$SAMPLE & done ``` Después de arreglar 1000 errores, en la carpeta AI4FOOD/READ_QC, se generan dos carpetas: V1_1065 y V1_1100, con los informes QC tanto antes como después de ejecutar FASTQC. **Para las dos muestras, ha tardado 1h:30m, luego parece tardar 45 min por muestra** :cold_sweat:. En conclusión, aunque tarde la misma vida para procesar cada muestra, **metaWRAP funciona** :smile_cat:, así que puedo ejecutarlo con el resto de muestras sin miedo. Paso las lecturas limpias de estas dos muestras a la carpeta **CLEAN_READS**: ```r=1 #dentro de la carpeta AI4FOOD mkdir CLEAN_READS for i in READ_QC/*; do b=${i#*/} mv ${i}/final_pure_reads_1.fastq CLEAN_READS/${b}_1.fastq mv ${i}/final_pure_reads_2.fastq CLEAN_READS/${b}_2.fastq done ``` Con este script obtengo cuatro ficheros: V1_1065_1.fastq, V1_1065_2.fastq, V1_1100_1.fastq y V1_1100_2.fastq, que se llaman igual que los archivos antes de ejecutar metaWRAP, pero son realmente los ficheros ya limpios. **Por eso, usarás los ficheros de CLEAN_READS para ejecutar MetaPhlAn4.** ::: Acortados los nombres de todos los archivos, he movido los ficheros de la carpeta "Lecturas_repetidas" a la carpeta "RAW_READS" usando FileZilla, de forma que al final tenemos en dicha carpeta dos ficheros .fq, con la misma estructura, para cada muestra. He cogido V1_1065 y V1_1100 de esa carpeta y los he movido a otra carpeta llamada "Reads ya hechos con metawrap", para que metaWRAP no me procese esas muestras que ya están procesadas en la carpeta "READ_QC". :::success **NÚMERO DE MUESTRAS: 95** (considerando 2141_a y 2141_b como dos muestras distintas). ::: **A fecha de 27/10/2022, he corrido metaWRAP con todas las muestras** (excepto las dos que ya he mencionado), es decir, **93 muestras (186 archivos .fq)**, usando: ```r=1 cd ~/victor/AI4FOOD for F in RAW_READS/*_1.fq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 1 -o READ_QC/$SAMPLE & done ``` :::danger **Lo dejaré corriendo todo el puente y a probar suerte** :crossed_fingers: ::: **A fecha del 6/11/2022, metaWRAP ha terminado con las 95 muestras** (aunque no conseguía hacerlas todas del tirón así que algunas tuve que hacerlas una a una). He movido todas las carpetas resultantes a la carpeta AI4FOOD/FINAL_READ_QC. Dentro de AI4FOOD, ejecuto finalmente el siguiente comando: ```r=1 #dentro de la carpeta AI4FOOD mkdir CLEAN_READS for i in FINAL_READ_QC/*; do b=${i#*/} mv ${i}/final_pure_reads_1.fastq CLEAN_READS/${b}_1.fastq mv ${i}/final_pure_reads_2.fastq CLEAN_READS/${b}_2.fastq done ``` Tras esto, dentro de la carpeta **AI4FOOD**, tenemos una carpeta llamada **CLEAN_READS** con las lecturas limpias de todas las muestras **(95 muestras - 190 archivos .fastq)**. ### **PASO 4: MetaPhlAn4** TASK: MetaPhlAn4 Ahora, vamos a ejecutar **MetaPhlan4** con todas nuestras muestras. Antes de hacerlo con todas ellas, vamos a hacer una prueba con un par: :::info **TEST ANTES DE USAR MetaPhlAn4 CON TODAS LAS MUESTRAS** Antes de ejecutar MetaPhlAn4 con todas las muestras, voy a hacerlo con un par de muestras (V1_1045 y V1_1065, que están en la carpeta "prueba_metaphlan"), para ver si todo va bien: ```r=1 conda activate mpa metaphlan --version #MetaPhlAn version 4.0.2 (22 Sep 2022) cd ~/victor/AI4FOOD/prueba_metaphlan touch metaphlan_script.pl #creo un script chmod a+x metaphlan_script.pl #para poder modificarlo vi metaphlan_script.pl #lo abro para modificarlo #copio lo siguiente en dicho fichero: #!/usr/bin/perl #dar de argumento solo los _1 foreach $ar (@ARGV) { # sacamos cada archivo: @field = split (/_/, $ar); $name = $ar; # name = nombre completo $name =~ s/_1.fastq//; # name = nombre antes de _ # obtenemos el _2 correspondiente: $r1 = $ar; # nombre completo de lectura con _1 $r2 = $ar; # lo mismo que r1 por ahora $r2 =~ s/1.fastq/2.fastq/; # cambia _1 por _2 # pasamos ambos a metaphlan: system("metaphlan $r1,$r2 --bowtie2db ~/victor/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt"); } ./metaphlan_script.pl *_1.fastq > metaphlanOUT_06_11_2022_17_50.txt 2>&1 #esto lo hago por si sale algún error. ``` **Todo salió correctamente aparentemente, así que lo voy a correr con el resto de muestras.** ::: Ejecutamos MetaPhlAn4 con todas las muestras restantes: ```r=1 conda activate mpa metaphlan --version #MetaPhlAn version 4.0.2 (22 Sep 2022) cd ~/victor/AI4FOOD/CLEAN_READS touch metaphlan_script.pl #creo un script chmod a+x metaphlan_script.pl #para poder modificarlo vi metaphlan_script.pl #lo abro para modificarlo #copio lo siguiente en dicho fichero: #!/usr/bin/perl #dar de argumento solo los _1 foreach $ar (@ARGV) { # sacamos cada archivo: @field = split (/_/, $ar); $name = $ar; # name = nombre completo $name =~ s/_1.fastq//; # name = nombre antes de _ # obtenemos el _2 correspondiente: $r1 = $ar; # nombre completo de lectura con _1 $r2 = $ar; # lo mismo que r1 por ahora $r2 =~ s/1.fastq/2.fastq/; # cambia _1 por _2 # pasamos ambos a metaphlan: system("metaphlan $r1,$r2 --bowtie2db ~/victor/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt"); } ./metaphlan_script.pl *_1.fastq > metaphlanOUT_06_11_2022_19_23.txt 2>&1 #esto lo hago por si sale algún error. ``` A fecha del 10/11/2022, MetaPhlan4 ha procesado correctamente todas las muestras, luego voy a fusionar los informes individuales de cada muestra en un único archivo: ```r=1 #Muevo los informes de las 2 muestras de prueba anteriores #a la carpeta CLEAN_READS, y entonces ejecuto el siguiente #comando. #En la carpeta victor/AI4FOOD/CLEAN_READS python ~/anaconda3/envs/mpa/bin/merge_metaphlan_tables.py profiled_*.txt > merged_abundance_table.txt ``` :::success **Ya tenemos la tabla de abundancias relativas generada por MetaPhlan4 para todas las muestra :cool:.** ::: ### **PASO 5: Análisis en R** TASK: En R, vamos a integrar los metadatos con las abundancias y vamos a sacar las diversidades :::warning **A fecha del 11/11/2022, vamos a usar como metadatos datos básicos (id, grupo, fecha nacimiento, peso, altura...). Más tarde iremos añadiendo más datos.** ::: **1) Vamos a generar la función metaphlanToPhyloseq, que nos será de utilidad para los análisis de diversidad.** ```r=1 metaphlanToPhyloseq <- function( tax, metadat=NULL, simplenames=TRUE, roundtointeger=FALSE, split="|"){ ## tax is a matrix or data.frame with the table of taxonomic abundances, rows are taxa, columns are samples ## metadat is an optional data.frame of specimen metadata, rows are samples, columns are variables ## if simplenames=TRUE, use only the most detailed level of taxa names in the final object ## if roundtointeger=TRUE, values will be rounded to the nearest integer xnames = rownames(tax) shortnames = gsub(paste0(".+\\", split), "", xnames) if(simplenames){ rownames(tax) = shortnames } if(roundtointeger){ tax = round(tax * 1e4) } x2 = strsplit(xnames, split=split, fixed=TRUE) taxmat = matrix(NA, ncol=max(sapply(x2, length)), nrow=length(x2)) colnames(taxmat) = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species", "Strain")[1:ncol(taxmat)] rownames(taxmat) = rownames(tax) for (i in 1:nrow(taxmat)){ taxmat[i, 1:length(x2[[i]])] <- x2[[i]] } taxmat = gsub("[a-z]__", "", taxmat) taxmat = phyloseq::tax_table(taxmat) otutab = phyloseq::otu_table(tax, taxa_are_rows=TRUE) if(is.null(metadat)){ res = phyloseq::phyloseq(taxmat, otutab) }else{ res = phyloseq::phyloseq(taxmat, otutab, phyloseq::sample_data(metadat)) } return(res) } ``` **2) Convertimos en objeto de phyloseq** Siguiendo como modelo el HackMD de Blanca y otro HackMD en el que se basa este a su vez, los pasos son los siguientes: * **[Metadatos](https://drive.google.com/file/d/14ncLjFVA_0ak10sAiadLdKyi7mOAlVo_/view?usp=sharing) que he importado** :::warning **De los metadatos que me pasó Blanca por correo, he eliminado: 2141 (es la muestra que está repetida, y como no sabía qué hacer con ella la he quitado en principio de los metadatos), 2144, 2093, 1257, 2134, 2158 y 2169 (que hubo problemas con ellas o no fueron secuenciadas).** ::: * **[Abundancias](https://drive.google.com/file/d/1NY6GVDobrLC_UWq1A0jVQnBAz6GR7wut/view?usp=sharing) que he importado** :::warning **De las abundancias he quitado 2141, ya que también la he quitado de los metadatos.** ::: ```r=1 #Importo los metadatos metadata_ai4food_v1<-read.csv('C:\\Users\\USUARIO\\Documents\\Archivos_imdea\\AI4FOOD\\metadatos_ai4food_v1.csv') #Indico los nombres de las filas row.names(metadata_ai4food_v1) <- metadata_ai4food_v1$id_voluntario sample.ai4food_v1 <- sample_data(metadata_ai4food_v1) #Como en MetaPhlan4 las muestras las tenemos como #profiled_V1_*, cambio el nombre de las filas sample_names(sample.ai4food_v1) <- gsub('^', 'profiled_V1_', + metadata_ai4food_v1$id_voluntario) #Importo las abundancias de MetaPhlan4 abundance_ai4food_v1<-read_csv("C:\\Users\\USUARIO\\Documents\\Archivos_imdea\\AI4FOOD\\abundancias_metaphlan_ai4food_v1.csv") rownames(abundance_ai4food_v1) <- abundance_ai4food_v1$clade_name abundances.ai4food_nocladename <- abundance_ai4food_v1 %>% select(-c(clade_name)) rownames(abundances.ai4food_nocladename)<-abundance_ai4food_v1$clade_name #Creo el objeto phyloseq phyloseq_ai4food_v1 = metaphlanToPhyloseq(abundances.ai4food_nocladename, metadat = sample.ai4food_v1) #hay que importar esta librería para usar rtree library(ape) #añadimos el árbol filogenético al objeto inicial random_tree=rtree(ntaxa(phyloseq_ai4food_v1), rooted = T, tip.label = taxa_names(phyloseq_ai4food_v1)) phyloseq_ai4food_v1.tree <- merge_phyloseq(phyloseq_ai4food_v1, random_tree) phyloseq_ai4food_v1.tree #phyloseq-class experiment-level object #otu_table() OTU Table: [ 3582 taxa and 93 samples ] #sample_data() Sample Data: [ 93 samples by 9 sample variables ] #tax_table() Taxonomy Table: [ 3582 taxa by 8 taxonomic ranks ] #phy_tree() Phylogenetic Tree: [ 3582 tips and 3581 internal nodes ] ``` :::success **¡Ya tenemos nuestro objeto de phyloseq!** :smile_cat: ::: --- **ANEXO: Evaluación de biomarcadores de diabetes y abundancias de los taxones identificados como importantes para PYM** Blanca me pasó una tabla de bioquímica para los voluntarios de AI4FOOD. He cogido las variables de HBA1C (en porc. y en mmol/mol), homocisteína, glucosa, HDL, PCR, HOMA y adiponectina, y puse esas columnas junto con la columna de individuos, creando el fichero que se importa a R a continuación: ```r=1 bioquimica_diabetes_ai4food<-read_csv("C:\\Users\\USUARIO\\Documents\\Archivos_imdea\\AI4FOOD\\bioquimica_biomarcadores_diabetes.csv") ``` **Antes de nada, vamos a aglomerar el objeto phyloseq de MetaPhlAn4 a nivel de género:** ```r=1 phyloseq_ai4food_v1_genus <- tax_glom(phyloseq_ai4food_v1.tree,taxrank="Genus",NArm=TRUE) ``` :::warning **OJO! El objeto phyloseq que se usa aquí de input es tras haber cambiado los metadatos que se indican a la función MetaphlanToPhyloseq(). Esto es porque he añadido a los metadatos las columnas de *bioquimica_diabetes_ai4food*, además de arreglar algunos valores (había columnas cuyo valor era 100 veces al real, como la hemoglobina glicosilada)**. ::: **A)** Ahora, vamos a ver cómo varían las abundancias en base al índice **HOMA-IR** de los participantes: ```r=1 tb_phyloseq_ai4food_genus<-psmelt(phyloseq_ai4food_v1_genus) #Lo que hago aquí abajo es porque #por algún motivo al transformar los datos de #metaphlan a phyloseq algunas abundancias se #multiplican por 1000, y eso pasa con la abundancia #de Dorea para el taxón 936. tb_phyloseq_ai4food_genus["5787","Abundance"]<-0.148923 #Divido en 4 cuartiles, siendo 1 las muestras #con las abundancias más bajitas y 4 las de las #abundancias más altas. tb_phyloseq_ai4food_genus_quantiles<-tb_phyloseq_ai4food_genus %>% mutate(quantile=ntile(homa,4)) tb_phyloseq_ai4food_genus_quantiles$quantile<-as.factor(tb_phyloseq_ai4food_genus_quantiles$quantile) #Me quedo solo con Phocea y Dorea tabla_taxones_importantes<-tb_phyloseq_ai4food_genus_quantiles[tb_phyloseq_ai4food_genus_quantiles$Genus %in% c("Phocea", "Dorea", "Lachnospira"), ] #Vemos las estadísticas para, por ejemplo, la HOMA-IR summary(tb_phyloseq_ai4food_genus_quantiles$homa) IQR(tb_phyloseq_ai4food_genus_quantiles$homa) ``` **Gráfica de los cuartiles para HOMA-IR y las abundancias de Phocea y Dorea (Murib. y P.NK3NB1_group no aparecen por diferencias entre SILVA y bowtie2):** ```r=1 ggplot(data = tabla_taxones_importantes, aes(x = quantile, y = Abundance, label = id_voluntario, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "red", show.legend = TRUE) + geom_text_repel( color = "white", # text color bg.color = "grey30", # shadow color bg.r = 0.15, # shadow radius )+ labs( title = "Abundances for Dorea, Lachnospira and Phocea", subtitle = "Groups made based on the HOMA-IR index", caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s",IQR(tb_phyloseq_ai4food_genus_quantiles$homa),quantile(tb_phyloseq_ai4food_genus_quantiles$homa, 0.25),quantile(tb_phyloseq_ai4food_genus_quantiles$homa, 0.75),quantile(tb_phyloseq_ai4food_genus_quantiles$homa, 0),quantile(tb_phyloseq_ai4food_genus_quantiles$homa, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5),plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black"), legend.background = element_rect(fill="lightblue",size=0.5, linetype="solid",colour ="darkblue")) ``` ![](https://i.imgur.com/VBpNDbu.png) **B) GLUCOSA** ```r=1 tb_phyloseq_ai4food_genus_quantiles_glucosa<-tb_phyloseq_ai4food_genus %>% mutate(quantile=ntile(glu_mg_dl,4)) tb_phyloseq_ai4food_genus_quantiles_glucosa$quantile<-as.factor(tb_phyloseq_ai4food_genus_quantiles_glucosa$quantile) tabla_taxones_importantes_glucosa<-tb_phyloseq_ai4food_genus_quantiles_glucosa[tb_phyloseq_ai4food_genus_quantiles_glucosa$Genus %in% c("Phocea", "Dorea", "Lachnospira"), ] ggplot(data = tabla_taxones_importantes_glucosa, aes(x = quantile, y = Abundance, label = id_voluntario, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "red", show.legend = TRUE) + geom_text_repel( color = "white", # text color bg.color = "grey30", # shadow color bg.r = 0.15, # shadow radius )+ labs( title = "Abundances for Dorea, Lachnospira and Phocea", subtitle = "Groups made based on glucose levels (mg/dl)", caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s",IQR(tb_phyloseq_ai4food_genus_quantiles_glucosa$glu_mg_dl),quantile(tb_phyloseq_ai4food_genus_quantiles_glucosa$glu_mg_dl, 0.25),quantile(tb_phyloseq_ai4food_genus_quantiles_glucosa$glu_mg_dl, 0.75),quantile(tb_phyloseq_ai4food_genus_quantiles_glucosa$glu_mg_dl, 0),quantile(tb_phyloseq_ai4food_genus_quantiles_glucosa$glu_mg_dl, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5),plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black"), legend.background = element_rect(fill="lightblue",size=0.5, linetype="solid",colour ="darkblue")) ``` ![](https://i.imgur.com/ZtesOau.png) **C) HBA1C (porcentaje)** ```r=1 tb_phyloseq_ai4food_genus_quantiles_hba1c<-tb_phyloseq_ai4food_genus %>% mutate(quantile=ntile(hba1c_perc,4)) tb_phyloseq_ai4food_genus_quantiles_hba1c$quantile<-as.factor(tb_phyloseq_ai4food_genus_quantiles_hba1c$quantile) tabla_taxones_importantes_hba1c<-tb_phyloseq_ai4food_genus_quantiles_hba1c[tb_phyloseq_ai4food_genus_quantiles_hba1c$Genus %in% c("Phocea", "Dorea", "Lachnospira"), ] ggplot(data = tabla_taxones_importantes_hba1c, aes(x = quantile, y = Abundance, label = id_voluntario, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "red", show.legend = TRUE) + geom_text_repel( color = "white", # text color bg.color = "grey30", # shadow color bg.r = 0.15, # shadow radius )+ labs( title = "Abundances for Dorea, Lachnospira and Phocea", subtitle = "Groups made based on HBA1C levels (percentage)", caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s",IQR(tb_phyloseq_ai4food_genus_quantiles_hba1c$hba1c_perc),quantile(tb_phyloseq_ai4food_genus_quantiles_hba1c$hba1c_perc, 0.25),quantile(tb_phyloseq_ai4food_genus_quantiles_hba1c$hba1c_perc, 0.75),quantile(tb_phyloseq_ai4food_genus_quantiles_hba1c$hba1c_perc, 0),quantile(tb_phyloseq_ai4food_genus_quantiles_hba1c$hba1c_perc, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5),plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black"), legend.background = element_rect(fill="lightblue",size=0.5, linetype="solid",colour ="darkblue")) ``` ![](https://i.imgur.com/QbQNtBz.png) **D) ADIPONECTINA** ```r=1 tb_phyloseq_ai4food_genus_quantiles_adip<-tb_phyloseq_ai4food_genus %>% mutate(quantile=ntile(adiponectina_ug_ml,4)) tb_phyloseq_ai4food_genus_quantiles_adip$quantile<-as.factor(tb_phyloseq_ai4food_genus_quantiles_adip$quantile) tabla_taxones_importantes_adip<-tb_phyloseq_ai4food_genus_quantiles_adip[tb_phyloseq_ai4food_genus_quantiles_adip$Genus %in% c("Phocea", "Dorea", "Lachnospira"), ] ggplot(data = tabla_taxones_importantes_adip, aes(x = quantile, y = Abundance, label = id_voluntario, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "red", show.legend = TRUE) + geom_text_repel( color = "white", # text color bg.color = "grey30", # shadow color bg.r = 0.15, # shadow radius )+ labs( title = "Abundances for Dorea, Lachnospira and Phocea", subtitle = "Groups made based on adiponectin levels (ug/ml)", caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s",IQR(tb_phyloseq_ai4food_genus_quantiles_adip$adiponectina_ug_ml),quantile(tb_phyloseq_ai4food_genus_quantiles_adip$adiponectina_ug_ml, 0.25),quantile(tb_phyloseq_ai4food_genus_quantiles_adip$adiponectina_ug_ml, 0.75),quantile(tb_phyloseq_ai4food_genus_quantiles_adip$adiponectina_ug_ml, 0),quantile(tb_phyloseq_ai4food_genus_quantiles_adip$adiponectina_ug_ml, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5),plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black"), legend.background = element_rect(fill="lightblue",size=0.5, linetype="solid",colour ="darkblue")) ``` ![](https://i.imgur.com/CWogX9G.png) **E) PROTEÍNA C REACTIVA** ```r=1 tb_phyloseq_ai4food_genus$pcr_mg_dl<-gsub("<","",tb_phyloseq_ai4food_genus$pcr_mg_dl) tb_phyloseq_ai4food_genus_quantiles_pcr<-tb_phyloseq_ai4food_genus %>% mutate(quantile=ntile(pcr_mg_dl,4)) tb_phyloseq_ai4food_genus_quantiles_pcr$quantile<-as.factor(tb_phyloseq_ai4food_genus_quantiles_pcr$quantile) tb_phyloseq_ai4food_genus_quantiles_pcr$pcr_mg_dl<-as.double(tb_phyloseq_ai4food_genus_quantiles_pcr$pcr_mg_dl) tabla_taxones_importantes_pcr<-tb_phyloseq_ai4food_genus_quantiles_pcr[tb_phyloseq_ai4food_genus_quantiles_pcr$Genus %in% c("Phocea", "Dorea", "Lachnospira"), ] ggplot(data = tabla_taxones_importantes_pcr, aes(x = quantile, y = Abundance, label = id_voluntario, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "red", show.legend = TRUE) + geom_text_repel( color = "white", # text color bg.color = "grey30", # shadow color bg.r = 0.15, # shadow radius )+ labs( title = "Abundances for Dorea, Lachnospira and Phocea", subtitle = "Groups made based on C reactive protein levels (mg/dl)", caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s",IQR(tb_phyloseq_ai4food_genus_quantiles_pcr$pcr_mg_dl),quantile(tb_phyloseq_ai4food_genus_quantiles_pcr$pcr_mg_dl, 0.25),quantile(tb_phyloseq_ai4food_genus_quantiles_pcr$pcr_mg_dl, 0.75),quantile(tb_phyloseq_ai4food_genus_quantiles_pcr$pcr_mg_dl, 0),quantile(tb_phyloseq_ai4food_genus_quantiles_pcr$pcr_mg_dl, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5),plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black"), legend.background = element_rect(fill="lightblue",size=0.5, linetype="solid",colour ="darkblue")) ``` ![](https://i.imgur.com/kIwKLUe.png) **F) IAS** ```r=1 tb_phyloseq_ai4food_genus$ias<-gsub("<","",tb_phyloseq_ai4food_genus$ias) tb_phyloseq_ai4food_genus_quantiles_ias<-tb_phyloseq_ai4food_genus %>% mutate(quantile=ntile(ias,4)) tb_phyloseq_ai4food_genus_quantiles_ias$quantile<-as.factor(tb_phyloseq_ai4food_genus_quantiles_ias$quantile) tb_phyloseq_ai4food_genus_quantiles_ias$ias<-as.double(tb_phyloseq_ai4food_genus_quantiles_ias$ias) tabla_taxones_importantes_ias<-tb_phyloseq_ai4food_genus_quantiles_ias[tb_phyloseq_ai4food_genus_quantiles_ias$Genus %in% c("Phocea", "Dorea", "Lachnospira"), ] ggplot(data = tabla_taxones_importantes_ias, aes(x = quantile, y = Abundance, label = id_voluntario, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "red", show.legend = TRUE) + geom_text_repel( color = "white", # text color bg.color = "grey30", # shadow color bg.r = 0.15, # shadow radius )+ labs( title = "Abundances for Dorea, Lachnospira and Phocea", subtitle = "Groups made based on C reactive protein levels (mg/dl)", caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s",IQR(tb_phyloseq_ai4food_genus_quantiles_ias$ias),quantile(tb_phyloseq_ai4food_genus_quantiles_ias$ias, 0.25),quantile(tb_phyloseq_ai4food_genus_quantiles_ias$ias, 0.75),quantile(tb_phyloseq_ai4food_genus_quantiles_ias$ias, 0),quantile(tb_phyloseq_ai4food_genus_quantiles_ias$ias, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5),plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black"), legend.background = element_rect(fill="lightblue",size=0.5, linetype="solid",colour ="darkblue")) ``` --- **3) Alfa-diversidad** :::info **12/11/22** Ahora que tenemos el objeto de phyloseq, vamos a hacer algunas gráfica de alfa-diversidad, por ver un poco cómo salen. **ESTAS NO SON LAS GRÁFICAS DEFINITIVAS**. ::: :::warning **MetaPhlan4 proporciona las abundancias relativas por millón, ¡luego hay que pasarlas a absolutas!** ::: ```r=1 GP1 <- transform_sample_counts(phyloseq_ai4food_v1.tree, function(x) 1E6 * x) ``` * **ALFA-DIVERSIDAD POR SEXO** :::info Con la función **geom_signif**, podemos hacer una Wilcoxon directamente para saber si hay diferencias significativas entre grupos. ::: :::info **Aquí no se muestra, pero geom_signif tiene una opción llamada map_signif_level = TRUE, por si en lugar de ver los p-valores queremos ver asteriscos.** ::: ```r=1 library(ggplot2) library(ggsignif) plot_richness(GP1,x="sexo",color="sexo", measures=c("Chao1","Shannon","Simpson"))+ geom_boxplot(aes(fill = sexo), alpha=.7) + theme_bw()+ geom_signif(comparisons = list(c("H", "M")), test = "wilcox.test") ``` ![](https://i.imgur.com/UcfHEcq.png) :::warning **En las mujeres hay algunas muestras con muy baja diversidad :shrug:, por lo demás no hay nada significativo.** ::: * **ALFA-DIVERSIDAD POR GRUPO** :::info **Recordemos que hay dos grupos. Los voluntarios del grupo 1 empezaron con la toma de datos "tradicional" (cuestionarios) y pasaron a la digital después, y los del grupo 2 lo hicieron al revés.** ::: ```r=1 #Para hacer las diversidades para los grupos, #he tenido que cambiar los nombres de los grupos #y llamarlos Grupo_1 y Grupo_2 sample_data(phyloseq_ai4food_v1.tree)$grupo<-gsub("^","Grupo_",sample_data(phyloseq_ai4food_v1.tree)$grupo) #aplico el cambio GP1 <- transform_sample_counts(phyloseq_ai4food_v1.tree, function(x) 1E6 * x) #Represento la gráfica plot_richness(GP1,x="grupo",color="grupo", measures=c("Chao1","Shannon","Simpson"))+ geom_boxplot(aes(fill = grupo), alpha=.7) + theme_bw()+ geom_signif(comparisons = list(c("Grupo_1", "Grupo_2")), test = "wilcox.test") ``` ![](https://i.imgur.com/eFOdqX6.png) :::warning **En ambos grupos, y según el índice, hay algunas muestras con muy baja diversidad :shrug:, por lo demás no hay nada significativo.** Vamos a ejecutar lo siguiente: ```r=1 alpha<-estimate_richness(GP1,measures=c("Chao1", "Shannon", "Simpson")) ``` Sale lo siguiente: ``` > alpha Chao1 se.chao1 Shannon Simpson profiled_V1_1045 458 0 1.1207841 0.66815498 profiled_V1_1065 593 0 3.2338593 0.92562897 profiled_V1_1077 881 0 4.2767708 0.96214035 profiled_V1_1093 551 0 2.8496152 0.86088154 profiled_V1_1100 679 0 2.6899304 0.80198145 profiled_V1_1101 620 0 2.1767978 0.80082966 profiled_V1_1109 534 0 4.2613279 0.96920066 profiled_V1_111 669 0 1.3902904 0.65033739 profiled_V1_1120 972 0 2.2447885 0.85737657 profiled_V1_1133 599 0 3.7191813 0.95079692 profiled_V1_1179 806 0 1.2870158 0.59800882 profiled_V1_1220 791 0 4.0930117 0.95531800 profiled_V1_1251 721 0 1.4642639 0.62432004 profiled_V1_1252 546 0 1.7907766 0.53316690 profiled_V1_1256 615 0 1.3384918 0.59836594 profiled_V1_1264 301 0 3.5804965 0.95033705 profiled_V1_1274 917 0 2.0769786 0.62450783 profiled_V1_1310 727 0 3.3493595 0.92135333 profiled_V1_1327 466 0 0.8054717 0.51720384 profiled_V1_1330 225 0 3.2766777 0.93575949 profiled_V1_1660 745 0 4.3973004 0.96749605 profiled_V1_1662 584 0 2.4635710 0.87929232 profiled_V1_1664 698 0 2.7456775 0.74983490 profiled_V1_1666 730 0 1.9782517 0.69641743 profiled_V1_1674 867 0 4.1162181 0.96574312 profiled_V1_1797 596 0 2.4283864 0.81546521 profiled_V1_1910 843 0 1.4731484 0.56965919 profiled_V1_1976 327 0 3.4798576 0.94633162 profiled_V1_1987 478 0 1.2697974 0.58442869 profiled_V1_1988 427 0 3.4835607 0.93806800 profiled_V1_1989 417 0 3.3832424 0.93630362 profiled_V1_1993 603 0 4.0085858 0.96147436 profiled_V1_2015 799 0 0.1629646 0.03451943 profiled_V1_2016 658 0 2.8003747 0.80498447 profiled_V1_2024 445 0 3.5969999 0.94984606 profiled_V1_2036 1015 0 2.2309315 0.85323910 profiled_V1_2081 770 0 2.2244384 0.79100473 profiled_V1_2119 614 0 3.4860098 0.94372631 profiled_V1_2130 89 0 3.3930027 0.94225657 profiled_V1_2131 469 0 2.5802092 0.77980024 profiled_V1_2132 449 0 2.7266042 0.85639840 profiled_V1_2133 845 0 2.7493030 0.80994723 profiled_V1_2136 303 0 3.3204443 0.92829385 profiled_V1_2137 740 0 1.9759039 0.76471393 profiled_V1_2138 802 0 3.9262640 0.95639276 profiled_V1_2139 311 0 3.5199108 0.94161193 profiled_V1_2140 356 0 3.5715895 0.94755219 profiled_V1_2142 563 0 1.1259524 0.66850573 profiled_V1_2143 812 0 3.8095518 0.95488416 profiled_V1_2145 380 0 2.2130936 0.80400708 profiled_V1_2146 331 0 2.8262230 0.82801482 profiled_V1_2148 786 0 4.0417166 0.95734153 profiled_V1_2149 637 0 1.8053571 0.67070653 profiled_V1_2150 752 0 2.0063301 0.79096060 profiled_V1_2151 465 0 1.9810015 0.77414176 profiled_V1_2152 656 0 2.8248762 0.86592518 profiled_V1_2153 764 0 3.9647738 0.95262790 profiled_V1_2154 698 0 2.3487403 0.65268042 profiled_V1_2155 766 0 3.8327683 0.94615068 profiled_V1_2156 455 0 3.1771693 0.91119814 profiled_V1_2157 757 0 3.9299199 0.95437163 profiled_V1_2159 603 0 3.1748375 0.88712981 profiled_V1_2160 391 0 3.3495494 0.93158599 profiled_V1_2161 476 0 2.2966779 0.84848940 profiled_V1_2162 338 0 3.5292608 0.94696325 profiled_V1_2163 575 0 3.6889053 0.95012541 profiled_V1_2164 404 0 3.1468619 0.90402906 profiled_V1_2165 889 0 1.8484146 0.74591232 profiled_V1_2166 707 0 0.9331350 0.43796136 profiled_V1_2167 505 0 3.6680590 0.95159036 profiled_V1_2168 317 0 3.5999058 0.94800132 profiled_V1_2170 628 0 0.9025977 0.52574977 profiled_V1_2171 664 0 4.0139274 0.95890517 profiled_V1_2172 841 0 4.3781260 0.96756931 profiled_V1_447 955 0 2.6299111 0.81640610 profiled_V1_564 360 0 1.1054131 0.31849786 profiled_V1_614 493 0 2.0991256 0.73674295 profiled_V1_644 944 0 1.6461217 0.63331626 profiled_V1_650 849 0 1.0338053 0.54509539 profiled_V1_664 813 0 3.0372094 0.87472456 profiled_V1_668 449 0 3.4708138 0.94036094 profiled_V1_752 839 0 4.0747638 0.95992392 profiled_V1_760 812 0 3.2905079 0.89213238 profiled_V1_770 387 0 3.0961896 0.92003015 profiled_V1_783 792 0 2.7504439 0.86448070 profiled_V1_787 567 0 0.7458913 0.50546845 profiled_V1_789 650 0 1.9161179 0.68500777 profiled_V1_814 762 0 4.2407292 0.96335188 profiled_V1_817 339 0 1.1179058 0.66801084 profiled_V1_819 636 0 3.6095459 0.94754351 profiled_V1_820 325 0 3.5863144 0.95033709 profiled_V1_922 500 0 1.0077230 0.34990259 profiled_V1_936 691 0 4.2242793 0.96809741 ``` **Según Chao1, hay dos muestras del Grupo 2 con una diversidad muy baja (por debajo de 250), que son las muestras 1330 y 2130. Sin embargo, esas mismas muestras, para Simpson y Shannon, tienen unos índices que están dentro de los boxplot, es decir, no son outliers**. **Para Shannon, hay una muestra del Grupo 1 con una diversidad muy baja (0.16), que es la muestra 2015. En este caso, el índice de Simpson coincide en este aspecto: también es muy baja la diversidad según este índice. Curiosamente, para Chao1 es 799, el cual no es un valor que lo identifique como outlier**. **Según Simpson, hay dos muestras del Grupo 2 con abundancias también bajas: 2166 y 564. En esto no coinciden los otros índices**. Como cada índice indica cosas distintas, ¿a lo mejor esto no tiene tanta relevancia? :shrug: ::: **4) 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) distances <- function(study_pseq) { for( i in dist_methods ){ # Calculate distance matrix iDist <- 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 <- adonis(dlist[[x]] ~ grupo, data = data.frame(sample_data(study_pseq))) return(as.data.frame(z$aov.tab)) }) 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") + geom_point(size = 2) + stat_ellipse(aes(group = grupo), linetype = 2) + theme_bw() + theme(plot.title = element_text(face = 'bold', size = 16), axis.title = element_text(size = 14), legend.title = element_text(size = 14), axis.text = element_text(size = 12), legend.text = element_text(size = 12)) 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 ``` ```r=1 library(cowplot) pseq.compositional_ai4food <- transform(phyloseq_ai4food_v1.tree, "compositional") ai4food_v1.beta <- distances(pseq.compositional_ai4food) ai4food_v1.beta.adonis <- adonis_calculator(ai4food_v1.beta$dlist, pseq.compositional_ai4food) ai4food_v1.beta.plot <- plotter_beta(dist_methods, pseq.compositional_ai4food, ai4food_v1.beta$pcoa_list) plist <- plotter_beta(dist_methods, pseq.compositional_ai4food, ai4food_v1.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') ``` **Resultado:** ![](https://i.imgur.com/wNnvuo9.png) :::info **LEYENDA** **a: unweighted UniFrac b: weighted UniFrac c: Bray-Curtis d: Jaccard** ::: **Resultados del PERMANOVA:** ``` ai4food_v1.beta.adonis $unifrac Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) grupo 1 0.1285365 0.1285365 1.17388 0.0127355 0.18 Residuals 91 9.9642388 0.1094971 NA 0.9872645 NA Total 92 10.0927754 NA NA 1.0000000 NA $wunifrac Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) grupo 1 0.2049563 0.2049563 0.9523937 0.01035746 0.495 Residuals 91 19.5833162 0.2152013 NA 0.98964254 NA Total 92 19.7882725 NA NA 1.00000000 NA $bray Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) grupo 1 0.3023115 0.3023115 0.8821877 0.009601291 0.631 Residuals 91 31.1842356 0.3426839 NA 0.990398709 NA Total 92 31.4865471 NA NA 1.000000000 NA $jaccard Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) grupo 1 0.3629074 0.3629074 0.9098754 0.009899648 0.628 Residuals 91 36.2957074 0.3988539 NA 0.990100352 NA Total 92 36.6586148 NA NA 1.000000000 NA ``` :::info **Si no entiendo mal, R^2 me indica el porcentaje de la varianza explicada por los grupos (aquí, es siempre <1% :sweat_smile:), y el p-valor si hay diferencias entre grupos (no hay nada significativo :sweat_smile:).** ::: ### **PASO 5A: Análisis en R (tras hacer la transformación SGB a GTDB)** Puesto que varios géneros, en MetaPhlan4, acaban llamándose GGB + 4 cifras, he hecho la transformación a GTDB para que sea más interpretable todo. Eso lo he hecho de esta manera: * Para pasar de lo que saca MetaPhlAn a la taxonomía GTDB (esto lo hice en el entorno mpa, y lo ejecuté en el directorio donde estén todos los "profiles" de metaphlan, creando primero el directorio gtdb): `for F in $(ls); do sgb_to_gtdb_profile.py -i $F -o gtdb/$F; done` * Para hacer el merge activé el entorno de conda metaphlan403 y me fui al directorio gtdb que acabamos de crear. Luego: `merge_metaphlan_tables.py profiled_* -o merged_abundance_table_gtdb.txt --gtdb_profiles` * Y al hacer el objeto phyloseq metí el parámetro simplenames = TRUE: `phyloseqin <- metaphlanToPhyloseq(data, metadata, simplenames = FALSE)` Partiendo del objeto *merged_abundance_table_gtdb.txt*, lo he transformado a objeto phyloseq, siguiendo las mismas instrucciones que puse más arriba. **Ya tenemos el objeto de phyloseq:** ``` > phyloseqin_ai4food_v1.tree phyloseq-class experiment-level object otu_table() OTU Table: [ 2343 taxa and 46 samples ] sample_data() Sample Data: [ 46 samples by 22 sample variables ] tax_table() Taxonomy Table: [ 2343 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 2343 tips and 2342 internal nodes ] ``` ```r=1 #Aglomeramos todo a nivel de género, eliminando además los NAs. phy_gen_ias <- tax_glom(phyloseqin_ai4food_v1.tree,taxrank="Genus",NArm=TRUE) #Saco las abundancias absolutas phy_gen_ias <- transform_sample_counts(phy_gen_ias, function(x) 1E6 * x) #Abundancias clr phyloseqin_ai4food_v1.tree_clr <- microbiome::transform(phy_gen_ias, "clr") #Transformación directamente a tabla de abundancias #relativas (1000 gracias a ChatGPT xddd) tb_clr <- psmelt(phyloseqin_ai4food_v1.tree_clr) tb_clr_relevant_stuff <- cbind(tb_clr$SampleID,tb_clr$Abundance,tb_clr$Genus) colnames(tb_clr_relevant_stuff) <- c("SampleID","Abundance","Genus") df_spread <- as_tibble(tb_clr_relevant_stuff) %>% spread(key = Genus, value = Abundance) ``` **Ahora, pasamos a calcular las diversidades (dividiendo a los participantes en dos grupos siguiendo el mismo criterio que para PYM y el IAS):** ```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_ias <- makeTreeSummarizedExperimentFromPhyloseq(phy_gen_ias) tse.list <- list("p_ias" = p_ias) 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_ias <- as.data.frame(colData(tse.list$PYM)) %>% select(chao1, shannon, simpson, pielou, relative, log_modulo_skewness, group_HEI) counter <<- 0 graphs_ias <-lapply(df_ias[ ,c("chao1", "shannon", "simpson", "pielou", "relative", "log_modulo_skewness")], function(a) {counter <<- counter + 1 ggplot(df_ias, aes(x = group_HEI, y = a)) + geom_boxplot(aes(fill = group_HEI), 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 = group_HEI), size = 1.5) + ylab(gsub('_', ' ', colnames(df_ias)[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_ias, common.legend = TRUE, legend = "right", align = "hv") ``` ![](https://i.imgur.com/I3naUqd.png) :::warning **Wilcoxon** para asegurarnos de que no hay diferencias significativas entre los dos grupos, antes y después de hacer la corrección FDR: ```r=1 wil.test <- bind_rows(pairwise_wilcox_test(df_ias, chao1 ~ group_HEI), pairwise_wilcox_test(df_ias, shannon ~ group_HEI), pairwise_wilcox_test(df_ias, simpson ~ group_HEI), pairwise_wilcox_test(df_ias, pielou ~ group_HEI), pairwise_wilcox_test(df_ias, relative ~ group_HEI), pairwise_wilcox_test(df_ias, log_modulo_skewness ~ group_HEI)) %>% adjust_pvalue(method = "BH") %>% add_significance() ``` ``` wil.test # 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 24 22 287 0.621 0.957 ns 2 shannon Group_1 Group_2 24 22 279 0.752 0.957 ns 3 simpson Group_1 Group_2 24 22 267 0.957 0.957 ns 4 pielou Group_1 Group_2 24 22 268 0.939 0.957 ns 5 relative Group_1 Group_2 24 22 248 0.736 0.957 ns 6 log_modulo_skewness Group_1 Group_2 24 22 324 0.193 0.957 ns ``` ::: Ahora, vamos a dibujar las **beta diversidades**: ```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]] ~ group_HEI, 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 = "group_HEI") + geom_point(size = 2) + stat_ellipse(aes(group = group_HEI), 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_ias_comp <- tax_glom(phyloseqin_ai4food_v1.tree,taxrank="Genus",NArm=TRUE) phy_gen_ias_comp.beta <- distances(phy_gen_ias_comp) phy_gen_ias_comp.beta.adonis <- adonis_calculator(phy_gen_ias_comp.beta$dlist, phy_gen_ias_comp) phy_gen_ias_comp.beta.plot <- plotter_beta(dist_methods, phy_gen_ias_comp, phy_gen_ias_comp.beta$pcoa_list) plist <- plotter_beta(dist_methods, phy_gen_ias_comp, phy_gen_ias_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') ``` ![](https://i.imgur.com/pPAyP5U.png) ``` > phy_gen_ias_comp.beta.adonis $unifrac Df SumOfSqs R2 F Pr(>F) group_HEI 1 0.1439022 0.02880172 1.304858 0.139 Residual 44 4.8524052 0.97119828 NA NA Total 45 4.9963074 1.00000000 NA NA $wunifrac Df SumOfSqs R2 F Pr(>F) group_HEI 1 0.1386128 0.02917669 1.322356 0.231 Residual 44 4.6121935 0.97082331 NA NA Total 45 4.7508063 1.00000000 NA NA $bray Df SumOfSqs R2 F Pr(>F) group_HEI 1 0.3370767 0.03864032 1.76851 0.067 Residual 44 8.3863664 0.96135968 NA NA Total 45 8.7234431 1.00000000 NA NA $jaccard Df SumOfSqs R2 F Pr(>F) group_HEI 1 0.4491828 0.03503853 1.597676 0.046 Residual 44 12.3704966 0.96496147 NA NA Total 45 12.8196794 1.00000000 NA NA ``` :::info **Llama la atención que, para Jaccard, el p-valor es significativo, pero el R^2 no es precisamente bueno xd** ::: --- **Tabla descriptiva de AI4FOOD** ```r=1 library(gtsummary) library(tidyverse) #¡OJO! Antes de hacer esto hay que cambiar los nombres #de las columnas (porque hay puntos feos) #usando colnames(datos_ias)[i] trial2 <- metadata_ai4food_v1 %>% select("Sex","HbA1c (%)","HbA1c_IFCC (mmol/mol)","Glucose (mg/dl)","C-reactive protein (mg/dl)",`HOMA-IR`,"Adiponectin (ug/ml)","group_HEI") %>% mutate(`HOMA-IR` = factor(`HOMA-IR`, levels = c("<1.96", "1.96 to 2.99", "\u22653"))) table2 <- tbl_summary( trial2, by = "group_HEI", # split table by group missing = "no", # don't list missing data separately label = "Adiponectin (ug/ml)" ~ "Adiponectin (\u03BCg/ml)" )%>% add_n() %>% # add column with total number of non-missing observations add_p(pvalue_fun = ~style_pvalue(.x, digits = 3))%>% # update the column header bold_labels() %>% add_q() %>% modify_header(label ~ "**Variable**") %>% modify_table_styling( columns = label, rows = label %in% "HbA1c (%)", footnote = "Legend: <5.7: Normal values 5.7 to 6.4: Prediabetes", text_format = "bold") %>% modify_table_styling( columns = label, rows = label %in% "HbA1c_IFCC (mmol/mol)", footnote = "Legend: \u226438: Normal values 39 to 47: Prediabetes", text_format = "bold") %>% modify_table_styling( columns = label, rows = label %in% "Glucose (mg/dl)", footnote = "Legend: \u226499: Normal values 100 to 125: Prediabetes", text_format = "bold") %>% modify_table_styling( columns = label, rows = label %in% "C-reactive protein (mg/dl)", footnote = "Legend: <1: Normal values \u22651: Abnormal values", text_format = "bold") %>% modify_table_styling( columns = label, rows = label %in% "HOMA-IR", footnote = "Legend: <1.96: Normal values 1.96 to 2.99: There might be insulin resistance \u22653: Insulin resistance", text_format = "bold") %>% modify_table_styling( columns = label, rows = label %in% "Adiponectin (\u03BCg/ml)", footnote = "Legend: <5: Unusual values \u22655: Abnormal values", text_format = "bold") %>% modify_caption("**Table 1. Summary statistics**") ``` ![](https://i.imgur.com/0HWw6cw.png) --- **Abundancias diferenciales AI4FOOD (solo para los individuos de los grupos):** ```r=1 #partimos del objeto phy_gen, con las abundancias absolutas #(son las únicas abundancias con las que esto funciona) phygen_prev <- subset_samples(phy_gen_ias) phygen_deseq <- phyloseq_to_deseq2(phygen_prev, ~ group_HEI) phygen_deseq <- DESeq(phygen_deseq, test="Wald", fitType="parametric") res_phygen <- results(phygen_deseq, cooksCutoff = FALSE, contrast=c("group_HEI", "Group_1", "Group_2")) #importante que esté la barra baja, si no no funciona alpha <- 0.05 coldtab <- res_phygen[which(res_phygen$pvalue < alpha), ] coldtab = cbind(as(coldtab, "data.frame"), as(tax_table(phygen_prev)[rownames(coldtab), ], "matrix")) #me quedo con aquellos cuyo log2FC sea mayor a 1 o menor a -1 subset_coldtab <- subset(coldtab, coldtab$log2FoldChange >=1 | coldtab$log2FoldChange <=-1) #quito los uncultured subset_coldtab <- subset_coldtab[!grepl("uncultured",subset_coldtab$Genus),] theme_set(theme_bw()) scale_fill_discrete <- function(palname = "Set1", ...) { scale_fill_brewer(palette = palname, ...) } # Genus order x = tapply(subset_coldtab$log2FoldChange, subset_coldtab$Genus, function(x) max(x)) x = sort(x, TRUE) subset_coldtab$Genus = factor(as.character(subset_coldtab$Genus), levels=names(x)) ggplot(subset_coldtab, aes(y=log2FoldChange, x=Genus, color = Phylum)) + geom_point(size=4) + theme(axis.text.x = element_text(angle=90, face = "bold.italic", color = "black", size = 11, hjust = 1, vjust = 0.5), axis.text.y = element_text(face="bold", color = "black"), axis.title.x = element_text(face="bold"), axis.title.y = element_text(face="bold"), legend.position = "right") + geom_hline(aes(yintercept=-1, linetype="dashed")) + geom_hline(aes(yintercept=1, linetype="dashed")) + guides(linetype = "none") + theme(legend.title = element_text(face = "bold"), legend.text = element_text(face = "bold.italic", size = 11)) + scale_y_continuous(n.breaks = 10) ``` :::warning **Esto es lo que se obtiene sin hacer más filtro aparte del del log2FC:** ![](https://i.imgur.com/xoO3X4O.png) **Esto es lo que se obtiene con pvalue < 0.05:** ![](https://i.imgur.com/hvdqMjI.png) **No hay ningún género con padj < 0.05 :sweat:** ::: ```r=1 compare_means(Abundance ~ quantile, data = important_taxa_table_hba1c[important_taxa_table_hba1c$Genus == "Bacteroides",]) ggplot(data = important_taxa_table_hba1c, aes(x = quantile, y = Abundance, label = SampleID, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "orange", show.legend = TRUE) + labs(caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s", IQR(tb_phyloseq_val_genus_quantiles_hba1c$hba1c_perc), quantile(tb_phyloseq_val_genus_quantiles_hba1c$hba1c_perc, 0.25), quantile(tb_phyloseq_val_genus_quantiles_hba1c$hba1c_perc, 0.75), quantile(tb_phyloseq_val_genus_quantiles_hba1c$hba1c_perc, 0), quantile(tb_phyloseq_val_genus_quantiles_hba1c$hba1c_perc, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5), plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold.italic", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black")) + stat_compare_means(comparisons = list(c(1,3)), label = "p.signif", size = 8, vjust = 0.5) ``` ```r=1 compare_means(Abundance ~ quantile, data = important_taxa_table_homa[important_taxa_table_homa$Genus == "Pseudoflavonifractor",]) ggplot(data = important_taxa_table_homa, aes(x = quantile, y = Abundance, label = SampleID, group = quantile, fill = quantile)) + geom_boxplot() + labs(x = "Quantile", y = "Abundance\n") + facet_wrap(~ Genus, scales = "free") + stat_summary(fun = "mean", geom = "point", aes(shape = "Mean"), size = 2, color = "orange", show.legend = TRUE) + labs(caption = sprintf("IQR: %s\n Q1: %s\t | Q3: %s\n Min: %s\t | Max: %s", IQR(tb_phyloseq_val_genus_quantiles_homa$homa.ir), quantile(tb_phyloseq_val_genus_quantiles_homa$homa.ir, 0.25), quantile(tb_phyloseq_val_genus_quantiles_homa$homa.ir, 0.75), quantile(tb_phyloseq_val_genus_quantiles_homa$homa.ir, 0), quantile(tb_phyloseq_val_genus_quantiles_homa$homa.ir, 1)), ) + theme( plot.title = element_text(color = "#0099f8", size = 18, face = "bold", hjust = 0.5), plot.subtitle = element_text(face = "bold.italic", hjust = 0.5), plot.caption = element_text(face = "bold.italic", size = 11, hjust = 1), strip.text = element_text(face="bold.italic", size=10), axis.title.x = element_text(size=12, face="bold", colour = "black"), axis.title.y = element_text(size=12, face="bold", colour = "black"), axis.text.x = element_text(size=10, face="bold", colour = "black"), axis.text.y = element_text(size=10, face="bold", colour = "black"), legend.title = element_text(size=11, face="bold", colour = "black")) + stat_compare_means(comparisons = list(c(1,3), c(3,4)), label = "p.signif", size = 8, vjust = 0.5) ``` ``` ```

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.