###### tags: `Estudiantes` `WGS` `Dieta`

# TFM Joaquín: "Gut microbiome as a marker for dietary intake in patients with celiac disease".

(i) Mining of public metagenomic datasets and dietary records.
(ii) Application of the Metagenomic Estimation of Dietary Intake (MEDI) tool to estimate dietary intake in a dataset of celiac disease patients (from the research group) and in publicly available control data.
(iii) Comparison of metagenomic dietary estimates with nutrient intake estimates derived from three-day dietary records.
(iv) Correlation of dietary intake data with metagenomic and metabolomic datasets.
(v) Implementation of machine learning approaches to identify dietary biomarkers by integrating metagenomic, metabolomic, and metagenomic-derived dietary intake estimates.
(vi) Biological interpretation of the results to elucidate their potential implications for public health.

## Temas de gestión:
- Ubicación en el Open Space planta 2 (Junto a Loles y Mar)
- Escribir correos (Con copia a mí y a Enrique) para:
- Pedir acceso a cluster UAM. :e-mail: [e-mail]
- Pedir acceso a cluster IMDEA. :e-mail: [e-mail]
- Pedir a Carlos TICS configuración de la VPN para poder conectarte desde tu casa. :e-mail: [e-mail]

:::warning
Si tienes alguna duda pregunta a Diego o Blanca. Pero en el correo sólo tienes que decir que eres estudiante del grupo y que necesitas acceso al cluster para hacer tu TFM
:::

## Organización tareas:

- Búsqueda bibliográfica de datasets públicos con datos de microbioma (WGS), datos de dieta (cuestionarios validados), metabolómica (opcional).
- La búsqueda bibliográfica debe seguir los pasos de PRISMA [ejemplo](https://journals.asm.org/doi/epub/10.1128/msystems.00143-25): 
- Buscar cualquier dataset público (pacientes sanos o enfermos, lo que haya) que tenga microbiota + dieta (preferiblemente recogida por cuestionarios validados).
- Puedes ver qué hay en el paquete de [CuratedMetagenomicData](https://www.bioconductor.org/packages/release/data/experiment/html/curatedMetagenomicData.html) por si nos sirve algo y luego descargar el raw data con su número de estudio en ENA o en donde esté publicado.
- Buscar en Pubmed, Scopus, etc.

- [Herramienta **MEDI**](https://pmc.ncbi.nlm.nih.gov/articles/PMC10871216/) - Hablar con Nico sobre la implementación y problemas que ha encontrado. Tu lo ejecutarías en los datos de celíacos (Datos de ECNR - `/shared/data/diego/kraken_celiac_noresp/data/WGS` | Datos públicos de controles `/shared/data/diego/kraken_celiac_control/data/tCD-TG`) y si encontramos otros datasets.

- **Datos de Dieta:** Hablar con Loles sobre los índices que ha calculado y preguntarle que necesitarías para implementarlo en el dataset de celíacos. Los datos que tenemos están [aquí](https://imdea-my.sharepoint.com/:f:/g/personal/judith_marcos_alimentacion_imdea_org/EuGGdZVv669NuNApLqrZaCABEZSZ-BuXczG_-5aNuF09Kg?e=TXiknu). Hay muchos archivos de datos, pero las muestras que necesitamos son las que están en el archivo `metadatos`. La idea sería poder trabajar con [este](https://jamesjiadazhan.github.io/dietaryindex_manual/index.html) paquete de R.

- No estaba en los objetivos primarios de tu TFM, pero es algo que tenemos pendiente de hacer... Qusiera que veas la herramienta [MetamobilePicker](https://metamobilepicker.readthedocs.io/en/master/) para caracterizar elemtntos móviles de resistencia dentro del metagenoma (lo podemos hablar a mi vuelta)- Aquí tendrás que hablar con Mar Morcillo para que te diga como está el estado de esto y como seguir. La idea es ejecutar este pipeline en el dataset de pacientes celíacos.

:::info
La primera semana tienes que centrarte en la búsqueda bibliográfica, lectura de artículos y depuración de datos de dieta, para tener claro lo que venga a futuro.
:::

# Trabajo con datos desacargados

- De todos los datasets: unificar metabolitos, metagenómica (Descargar de `curatedmetagenomicdata` ver tutorial) valorar correción de batch effect con Mmuphin, datos dieta.
- Descargar datos crudos para hacer MEDI. Los datos crudos de celiacos te los paso yo.
- Ver índices que ha calculado Loles (Te dí acceso a la carpeta).
- Primero hacer análisis exploratorio de todos los datos.
- Ejecutar MEDI
- Pensar de qué manera podemos comprar MEDI con cuestionarios (índices que calule etc.)

# Metamobilepicker

:warning: Hacer en el cluster de la UAM.

- Ejecutar los trabajos en: `/home/proyectos/imdeaalim/jgcorder`
- Crear un directorio con las muestras que vas a usar.
- Crear archivo config primero, que debe ser algo así: (Este es para una sóla muestra, la idea es que puedas mandar varias a la vez, aquí tienes que ver como se hace)

datadir: /lustre/NodoBIO/imdeaalim/ljmarcos
host: /lustre/local/metamobilepicker/0.7.2/lib/python3.10/site-packages/MetaMobilePicker/test/no_h$
max_mem: 8
outdir: /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta
samples:
  V1_1077:
    fwd: /lustre/NodoBIO/imdeaalim/ljmarcos/clean_data/V1_1077_1.fastq
    rev: /lustre/NodoBIO/imdeaalim/ljmarcos/clean_data/V1_1077_2.fastq
threads:
  big: 8
  huge: 16
  medium: 8
  small: 4

- Ejecutar metamobile usando el script de slurm similar a esto, Ejemplo `script.sh` modifica tu email, los directorios de output y error que es donde se guardan los logs.

#!/bin/bash
#SBATCH -p bioinfo                      #partition/queue name
#SBATCH --job-name=metamobile1                #Job name
#SBATCH -N 1                   #Un nodo
#SBATCH -n 16                   #Lo que tiene por default metamobile
#SBATCH -t 20:00:00             #Time limit hrs:min:sec
#SBATCH --output=/home/ljmarcos/metamobile/log-%j.o #Log de salida
#SBATCH --error=/home/ljmarcos/metamobile/log-%j.e #Log de errores
#SBATCH --mail-user=[e-mail]
#SBATCH --mail-type=END,FAIL

# Copiar a temporal los ficheros de entrada necesarios para metamobile
mkdir  /temporal/$SLURM_JOB_USER/$SLURM_JOB_ID
cp -r data /temporal/$SLURM_JOB_USER/$SLURM_JOB_ID
cd   /temporal/$SLURM_JOB_USER/$SLURM_JOB_ID
mkdir  output

# Ejecutar el comando Metamobile
module load metamobilepicker/0.7.2
metamobilepicker run -c /home/proyectos/imdeaalim/ljmarcos/V1_1077_config.yaml

# Mover a imdeaalim los resultados del cálculo.
mv output /lustre/NodoBIO/imdeaalim/ljmarcos/output

- Lanzar el script de SLURM.

sbatch -A imdeaalim_serv -p bioinfo script.sh 

:memo: La última prueba que he hecho me dió este error:

    # Filters contigs based on length, with a minimum of either 1000 or 500 bp.
        echo "Filtering assembly on 1000bp" > /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/log/assembly/V1_1077_length_filtering_out.txt
        seqkit seq -m 1000 /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/V1_1077_metaspades/contigs.fasta > /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/V1_1077_metaspades/V1_1077_1000bp.fasta 2> /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/log/assembly/V1_1077_length_filtering_error.txt
        
        echo "Moving assembly to temp directory for plasclass" >> /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/log/assembly/V1_1077_length_filtering_out.txt
        mkdir -p tmp 
        cp /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/V1_1077_metaspades/V1_1077_1000bp.fasta tmp/V1_1077_1000bp.fasta 2>> /lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/log/assembly/V1_1077_length_filtering_error.txt

        
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Removing output files of failed job cut_assembly since they might be corrupted:
/lustre/NodoBIO/imdeaalim/ljmarcos/out_meta/V1_1077/V1_1077_metaspades/V1_1077_1000bp.fasta
[Wed Jul  2 14:36:10 2025]
Finished job 12.
5 of 21 steps (24%) done
[Wed Jul  2 14:43:20 2025]
Finished job 14.
6 of 21 steps (29%) done
Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
[2025-07-02 14:43 CRITICAL] Command 'snakemake --use-conda --snakefile /lustre/local/metamobilepicker/0.7.2/lib/python3.10/site-packages/MetaMobilePicker/Snakefile --cores 16         --configfile /home/proyectos/imdeaalim/ljmarcos/V1_1077_config.yaml --rerun-incomplete  --nolock --latency-wait 60 --directory /lustre/local/metamobilepicker/0.7.2/lib/python3.10/site-packages/MetaMobilePicker --use-singularity        --singularity-prefix /lustre/local/metamobilepicker/0.7.2/lib/python3.10/site-packages/MetaMobilePicker/.snakemake --singularity-args '--home $PWD/.. --bind $TMPDIR,$PWD/..'' returned non-zero exit status 13.

# Phyloseq

:::danger
Lo ideal es que juntes el objeto phyloseq que te pasé con los datos de curated que has seleccionado y ejecutar Mmuphin para corregir batch effect. OJO si usas counts o datos transformados.
:::

## Análisis de diversidades
:::info
El análisis de diversidad puedes hacerlo con `phyloseq` o con `microbiome`
:::

De aquí sale el objeto phyloseq que te pasé:

#Loading libraries
library(tidyverse)
library(phyloseq)
library(curatedMetagenomicData)
library(ape)
library(microbiome)
library(factoextra)

#Abundance table
Abundance <- read_csv("Data/Metagenomica/merged_abundance_table.txt")
Abundance <- column_to_rownames(Abundance, var="clade_name")
Abundance$NCBI_tax_id <- NULL

#Converte to phyloseq object
#Metadata after clustering to add cluster as a variable
DATA_plus <- read_csv("Data/Objetos/data_plus.csv")
DATAmfa.cluster <- read_csv("Data/Objetos/DATAmfa.cluster.csv")
DATA <- DATAmfa.cluster %>% select(ID,Cluster)
names(DATA) <- c("Sample_ID", "Cluster")
DATA <- DATA %>% left_join(DATA_plus)
DATA$Cluster <- factor(DATA$Cluster, levels = c("1", "2"), labels=c("1","2"))
row.names(DATA) <- DATA$Sample_ID
sample <- sample_data(DATA)
sample_names(sample) <- DATA$Sample_ID

#Load script Function metaphlantophyloseq
phyloseqin= metaphlanToPhyloseq(Abundance, metadat = sample)
phyloseqin

#Tree 
random_tree=rtree(ntaxa(phyloseqin), rooted = T, tip.label = taxa_names(phyloseqin))
#Objeto phyloseq con arbol
physeq.tree <- merge_phyloseq(phyloseqin, random_tree)

### Alfa diversidad
En los datos generados con metaphlan :custard:
Necesitamos transformar a conteos absolutos para las diversidades. Porque metaphlan sólo dá abundancias relatvas por millon.
Para eso multiplicamos la abundancia de cada taxón por un millón - conteos por millón.
Debido a cómo funciona MetaPhlAn, no es posible tener "read counts". Sin embargo, se pueden tener "pseudo" conteos multiplicando las abundancias relativas por una constante y redondeando al número entero más cercano.

GP1 <- transform_sample_counts(physeq.tree, function(x) 1E6 * x)

En nuestro caso:

#Transfom counts for diversity analysis
NRCD <- transform_sample_counts(physeq.tree, function(x) 1E6 * x)
plot_richness(NRCD)
rich <- richness(NRCD)
plot_rich <- plot_richness(NRCD, color = "Cluster", x = "Cluster", measures = c("Chao1", "Simpson", "Shannon")) + geom_boxplot(aes(fill = Cluster), alpha=.7) + scale_color_manual(values = c("#B2ABD2", "#b2df8a")) + scale_fill_manual(values = c("#B2ABD2", "#b2df8a")) + theme_bw()

### Beta diversidad

#Compositional analisis transform to CLR 
pseq.compositional <- transform(physeq.tree, "compositional")

#Unifrac
ordinated_taxa_unifrac = ordinate(physeq.tree, method="MDS", distance="unifrac",weighted=TRUE)
plot_ordination(physeq.tree, ordinated_taxa_unifrac, color="Cluster", title = "Unifrac MDS Analysis") + theme_bw()

### Abundancias

Aquí ya ve [este](https://microbiome.github.io/OMA/docs/devel/pages/beta_diversity.html) tutorial y escoge los análisis que te interesen más.