###### tags: `16S rRNA` `WGS` `tortugas` `tutorials` # Para Helena :turtle: [ToC] ## Tutoriales y webs utiles * Tutorial qiime2 - Curso EMBO-Metagenomics. (Giuseppe D'Auria ) [Descargar](https://drive.google.com/file/d/1zpgpkPF7PtdLu41iLcbWO0flbJrK4QJo/view?usp=sharing) * Tutorial comandos linux. (Giuseppe D'Auria ) [Descargar Intro](https://drive.google.com/file/d/12bwKDoSDHbUDMgBspDnSvcX1auZMGGgC/view?usp=sharing). [Descargar tutorial](https://drive.google.com/file/d/18g4uhlencY7GlcrxrYWtclnjK-z6BD2h/view?usp=sharing). * [Biobakery ](https://github.com/biobakery/biobakery)tools :cake: * [Phyloseq](https://joey711.github.io/phyloseq/) * [Microbiome package](https://microbiome.github.io/tutorials/). ## Instalaciones - Instalar[METAWRAP](https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md) :memo: Usar lo siguiente como guia, modifica lo necesario, ojo con la versiones y actualizaciones. **Mejor simpre lo más actualizado.** Mejor vé [aquí](https://github.com/bxlab/metaWRAP) la instalación de Metawrap, que está más completo. - Instalar [Metaphlan](https://github.com/biobakery/MetaPhlAn/wiki/MetaPhlAn-4.1) - Instalar [Kraken](https://github.com/DerrickWood/kraken2/wiki/Manual) ## Conectarte al cluster En el terminal escribes: `ssh usuario@login2.ccc.uam.es ` Para pasar los datos a tu carpeta de ejecución: `/home/proyectos/imdeaalim/compubio/helena` En donde estan los datos de tus tortuguitas escribes: `scp *.fastq usuario@/home/proyectos/imdeaalim/compubio/helena/WGS` ## Análisis metagenomas WGS --- <font color = 'steelblue'> <p style="text-align:center;"><b><i>Comentarios de Blanca</i></b></p> </font> #### Comprobación de los ficheros Para asegurarnos de que la descarga ha ido bien, comparamos los **códigos md5** de los ficheros descargados con los que venían en el ENA. Para ello: - Sacamos los archivos md5 de cada fastq con este comando (por ahora **sigo con los `.gz`**): ```bash for file in $(find -mindepth 2 -type f); do echo $file; md5sum $file > $file.md5; done # pongo el echo para poder seguir el proceso ``` - Lo concatenamos todo en un solo fichero: ```bash for file in $(find -mindepth 2 -type f | grep md5); do cat $file >> md5_after_unzipping.txt ; done ``` - Ahora tenemos que generar el fichero `md5_paths_from_ENA.txt`, que tiene esta pinta: ``` 9676cfdd5e318bcbe70219d314089925 ./ERR414230/ERR414230_1.fastq.gz 3ee413f2e6cb17dbf5deac241e19f512 ./ERR414231/ERR414231_1.fastq.gz be2a52a759d8cd172714494d1eeb8f69 ./ERR414232/ERR414232_1.fastq.gz df45c11ad71f7b5c1ec5d7f6eac08ff2 ./ERR414233/ERR414233_1.fastq.gz c9d2db34f482833801f9df5535b3be86 ./ERR414234/ERR414234_1.fastq.gz b27940088b9b8df21396f39d986f3637 ./ERR414235/ERR414235_1.fastq.gz ``` > - En este caso, yo generé el fichero así: > ```bash >cat subsetENA_ftp_md5.tsv | awk '{print $4 " " $5}' > tmp45.txt #este archivo tiene solamente los archivos que queremos >wc -l tmp45.txt >tail -n 176 tmp45.txt > tmp_nohead.txt # queremos evitar encabezado >rm tmp45.txt >sed 's/ftp.sra.ebi.ac.uk\/vol1\/fastq\/ERR..././g' tmp_nohead.txt > tmp_paths.txt # elimina parte del path que no queremos >rm tmp_nohead.txt >grep -v ';' tmp_paths > part1.txt >grep ';' tmp_paths | tr ' ' ';' > semicolon.txt # columnas 1,3 son read1 y 2,4 son read2 >cut -d ';' -f1,3 semicolon.txt | tr ';' ' ' > read1.txt >cut -d ';' -f2,4 semicolon.txt | tr ';' ' ' > read2.txt >cat read1.txt read2.txt part.txt > md5_paths_from_ENA.txt >cp md5_paths_from_ENA.txt samples >cd samples >``` > Para las muestras del PRJEB2054: > ```bash >cat filereport_read_run_PRJEB2054_tsv.txt | grep 'bgi-MH' | cut -f29 > fastq_md5.txt #obtiene columna fastq_md5 >cat filereport_read_run_PRJEB2054_tsv.txt | grep 'bgi-MH' | cut -f30 > fastq_ftp.txt #columna fastq_ftp >paste fastq_md5.txt fastq_ftp.txt > both.txt #aqui podemos eliminar los dos .txt anteriores >sed -i 's/ftp.sra.ebi.ac.uk\/vol1\/fastq\/ERR..././g' both.txt >cut -d ';' -f1,4 both.txt | tr ';' ' ' > read.txt >cut -d ';' -f2,5 both.txt | tr ';' ' ' > read1.txt >cut -d ';' -f3,6 both.txt | tr ';' ' ' > read2.txt >cat read.txt read1.txt read2.txt > md5_paths_from_ENA.txt >mv md5_paths_from_ENA.txt RAW_READS >cd RAW_READS >``` :::success Finalmente, el comando que corremos es: ```bash md5sum --check md5_paths_from_ENA.txt > md5sum_check_log.txt 2>&1 ``` En el fichero `md5sum_check_log.txt` podremos ver si los ficheros se descargaron bien: ``` ./ERR414230/ERR414230_1.fastq.gz: OK ``` O si no fue el caso: ``` ./ERR414241/ERR414241_1.fastq.gz: FAILED ``` ::: > Para las muestras del PRJEB5024 me salen varios `FAILED open or read` porque me olvidé de **excluir las muestras con `lane`** de los MH0006 y MH0012. Con una comprobación rápida podemos ver que los FAILED coinciden con estos ficheros y que no hay de qué preocuparse: > ```bash > cat filereport_read_run_PRJEB2054_tsv.txt | grep 'lane' | cut -f6 #vemos los ERR de esos ficheros > cat md5sum_check_log.txt | grep 'FAILED open' | cut -d '/' -f2 | sort | uniq # OK! > ``` > Ha habido tres que sí que han fallado: > ```bash > cat md5sum_check_log.txt | grep 'FAILED' | grep -v open > ./ERR011189/ERR011189_1.fastq.gz: FAILED > ./ERR011088/ERR011088_2.fastq.gz: FAILED > ./ERR011229/ERR011229_2.fastq.gz: FAILED > ``` > Saco las ftp para intentar descargarlos con `wget`: > - ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011189/ERR011189_1.fastq.gz > - ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011088/ERR011088_2.fastq.gz > - ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011229/ERR011229_2.fastq.gz > > Tampoco ha funcionado, revisar qué pasa aquí. De momento las guardo en un directorio aparte y ya pensaré en ellas. > >Por otro lado, en estas muestras se da el problema de que tenemos un tercer fichero con el que no sé qué pasa. > - Mirando con el run ERR011087, los \_1 y \_2 tienen 46560224 líneas, mientras que el tercero en discordia tiene 3680 (12652 veces menos líneas). Si hago grep con el id del run, son 11640056 lecturas frente a 920 (lo mismo entre 4). > - Mirando los ID, haciendo head, tail, etc., no veo que los reads del tercer fichero aparezcan en los paired. Creo que son lecturas que no parearon con los otros dos ficheros, ahora queda **decidir qué hacemos con ellas.** > - De momento voy a correr metaWRAP con las paired y ya veremos qué pasa con las otras. #### Configuración del entorno e instalaciones Para la instalación seguí las [instrucciones del github](https://github.com/bxlab/metaWRAP#installation) Una cosa que tuve que hacer durante la instalación fue **descargar e indexar el hg38**. Seguí las [instrucciones de metaWRAP](https://github.com/bxlab/metaWRAP/blob/master/installation/database_installation.md#making-host-genome-index-for-bmtagger). El path que puse a la DB en el fichero config-metawrap es el siguiente: ``` BMTAGGER_DB=/home/lmarcos/BMTAGGER_INDEX ``` Los pasos que seguí para crear el entorno fueron: ```bash=v conda create --name metawrap-2 python=2.7 conda activate metawrap-2 conda config --add channels defaults conda config --add channels conda-forge conda config --add channels bioconda conda config --add channels ursky conda install --only-deps -c ursky metawrap-mg # tarda un poco pero lo saca ``` <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fin comentarios de Blanca</i></b></p> </font> ------- ## Empezar análisis con KRAKEN2 1. Construimos la libreria usando todas las secuencias disponibles en NCBI de virus, bacterias, hongos y protozoos ```bash = 1 #!/bin/bash #SBATCH -p bioinfo #partition/queue name #SBATCH --job-name=dbtest #Job name #SBATCH -N 1 #Un nodo #SBATCH -n 15 #15 cores #SBATCH -t 08:00:00 #Time limit hrs:min:sec #SBATCH --output=/home/ljmarcos/kraken/log-%j.o #Log de salida #SBATCH --error=/home/ljmarcos/kraken/log-%j.e #Log de errores #SBATCH --mail-user=judith.marcos@alimentacion.imdea.org #SBATCH --mail-type=ALL # creating variables for folders KRAKEN_DB=(/home/proyectos/imdeaalim/ljmarcos/kraken/kraken_db) ## Go to database folder cd $KRAKEN_DB #Cargar el modulo module load kraken2/2.1.3 ## Download taxonomy kraken-build --download-taxonomy --db kraken_db ## Download and install standard genomic libraries echo Downloading standard Kraken libraries kraken-build --download-library bacteria --db kraken_db kraken-build --download-library viral --db kraken_db kraken-build --download-library archaea --db kraken_db kraken-build --download-library fungi --db kraken_db kraken-build --download-library protozoa --db kraken_db echo Done with downloading and installing standard Kraken libraries ## Build the actual database (from libraries) kraken-build --build --db kraken_db --threads 15 ## Remove all intermediate files kraken-build --clean --db kraken_db --threads 15 echo Checking files ls echo Finished ``` 2. Ejecutamos KRAKEN2 ```bash =1 # Define los archivos de salida kraken_output="${kraken_output_dir}/${sample}_kraken_output.txt" kraken_report="${kraken_report_dir}/${sample}_kraken_report.txt" # Ejecuta Kraken2 kraken2 --db res/fungi_db --paired --confidence 0.1 --output $kraken_output --report $kraken_report $read1 $read2 ``` 3. Resultados con KRAKEN2 Primero pasamos los archivos de kraken a formato biom para que podamos importarlos a phyloseq. SE neecsita el programa [kraken-biom](https://github.com/smdabdoub/kraken-biom). ```bash kraken-biom *.txt --fmt json -o Kraken_completo.biom ``` 4. Análisis en R ```r=1 library(phyloseq) library(vegan) library(ggplot2) library(dplyr) library(RColorBrewer) library(patchwork) library(tidyr) library(tibble) library(tidyverse) library(ape) library(microbiome) library(genefilter) #Importar datos data <- import_biom("Data/Kraken_completo.biom") #Añadir metadatos metadata <- read_csv("Data/metadata.csv") row.names(metadata) <- metadata$SampleID sample <- sample_data(metadata) sample_names(sample) <- metadata$SampleID #Limpieza de datos data@tax_table@.Data <- substring(data@tax_table@.Data, 4) colnames(data@tax_table@.Data) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species") #Creación del árbol random_tree <- rtree(ntaxa(data), rooted = T, tip.label = taxa_names(data)) #Añadir arbol data <- merge_phyloseq(data, random_tree, sample) #Aglomerar a nivel de especie data <- tax_glom(data, "Species") #Filtrar #Hay que filtrar porque aquellas muestras con contajes muy bajos pueden ser artefactos de la secuenciación #Compositional analísis transform to CLR pseq.compositional <- transform(data, "compositional") #Escogemos como filtros aquellos taxones con una abundanica >5E-6 en al menos una muestra para arqueas, eucariotas y viruses # k = 1; A = 5E-6 flist <- filterfun(kOverA(1, 5E-6)) #Escogemos como filtros aquellos taxones con una abundanica >5E-3 en al menos una muestra para bacterias # k = 1; A = 5E-3 flist2 <- filterfun(kOverA(1, 5E-3)) # Filter taxa physeq.fil <- filter_taxa(pseq.compositional, flist) physeq.fil2 <- filter_taxa(pseq.compositional, flist2) physeq.fil <- prune_taxa(physeq.fil, pseq.compositional) physeq.fil.bacteria <- prune_taxa(physeq.fil2, pseq.compositional) #Plot de todos los reinos de bichos encontrados library(randomcoloR) #Preparar datos tb_sp <- psmelt(physeq.fil) #Crear variable con nombre completo de la especie tb_name <- tb_sp %>% unite(Genus,Species,col="Species_complete",sep=" ") tb_sp <- tb_sp %>% left_join(tb_name) #COn bacterias porque filtramos más tb_sp.bac <- psmelt(physeq.fil.bacteria) #Crear variable con nombre completo de la especie tb_name.bac <- tb_sp.bac %>% unite(Genus,Species,col="Species_complete",sep=" ") tb_sp.bac <- tb_sp.bac %>% left_join(tb_name.bac) #Virus tb_virus <- tb_sp %>% filter(Kingdom == "Viruses") paleta42 <- distinctColorPalette(k = 42, altCol = FALSE, runTsne = FALSE) plot_virus <- ggplot(tb_virus, aes(Sample, Abundance ,fill=Species_complete)) + geom_col(position="fill") + scale_fill_manual(values=paleta42) Plot_virus <- plot_virus + theme_minimal() + xlab("Sample") + 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")) Plot_virus + theme(legend.position = "bottom") + coord_flip() #Archaea tb_archaea <- tb_sp %>% filter(Kingdom == "Archaea") paleta22 <- distinctColorPalette(k = 22, altCol = FALSE, runTsne = FALSE) plot_archaea <- ggplot(tb_archaea, aes(Sample, Abundance ,fill=Species_complete)) + geom_col(position="fill") + scale_fill_manual(values=paleta22) Plot_archaea <- plot_archaea + theme_minimal() + xlab("Sample") + 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")) Plot_archaea + theme(legend.position = "bottom") + coord_flip() #Eukaryota #Diferenciar hongos y protistas tb_eukaryota <- tb_sp %>% filter(Kingdom == "Eukaryota") #Hongos tb_fungi <- tb_eukaryota %>% filter(Phylum == "Ascomycota" | Phylum == "Basidiomycota" | Phylum == "Microsporidia" | Phylum == "Mucormycota") paleta44 <- distinctColorPalette(k = 44, altCol = FALSE, runTsne = FALSE) plot_fungi<- ggplot(tb_fungi, aes(Sample, Abundance ,fill=Species_complete)) + geom_col(position="fill") + scale_fill_manual(values=paleta44) plot_fungi <- plot_fungi + theme_minimal() + xlab("Sample") + 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")) plot_fungi + theme(legend.position = "bottom") + coord_flip() #Diferenciar hongos y protistas tb_eukaryota <- tb_sp %>% filter(Kingdom == "Eukaryota") #Protistas tb_protistas <- tb_eukaryota %>% filter(Phylum == "Apicomplexa" | Phylum == "Bacillariophyta" | Phylum == "Euglenezoa" | Phylum == "Evosea" | Phylum == "Fornicata" | Phylum == "Oomycota" | Phylum == "Parabasalia") paleta18 <- distinctColorPalette(k = 18, altCol = FALSE, runTsne = FALSE) plot_protistas <- ggplot(tb_protistas, aes(Sample, Abundance ,fill=Species_complete)) + geom_col(position="fill") + scale_fill_manual(values=paleta18) plot_protistas <- plot_protistas + theme_minimal() + xlab("Sample") + 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")) plot_protistas + theme(legend.position = "bottom") + coord_flip() #Bacteria tb_bacteria <- tb_sp.bac %>% filter(Kingdom == "Bacteria") paleta100 <- distinctColorPalette(k = 100, altCol = FALSE, runTsne = FALSE) plot_bacteria<- ggplot(tb_bacteria, aes(Sample, Abundance ,fill=Species_complete)) + geom_col(position="fill") + scale_fill_manual(values=paleta100) plot_bacteria <- plot_bacteria + theme_minimal() + xlab("Sample") + 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")) plot_bacteria + theme(legend.position = "bottom") + coord_flip() #Bacteria por género: plot_bacteriaG <- ggplot(tb_bacteria, aes(Sample, Abundance ,fill=Genus)) + geom_col(position="fill") + scale_fill_manual(values=paleta100) plot_bacteriaG <- plot_bacteriaG + theme_minimal() + xlab("Sample") + 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")) plot_bacteriaG + theme(legend.position = "bottom") + coord_flip() ``` ## Análisis con [metaphlan4](https://github.com/biobakery/biobakery/wiki/metaphlan4) 1. Instalar metaphlan4 ` conda install -c bioconda metaphlan` 2. Ejecutar Metaphlan4 # Análisis en R :computer: :::info `Phyloseq` es un paquete de R que se usa mucho para análisis de microbioma, los datos de curated se pueden pasar fácilmente. Un objeto `Phyloseq` está formado por tres elementos necesarios: **1. La tabla de taxonomía. 2. La tabla de metadata. 3. La tabla de abundancia.** Y además... 4. La información filogenética. No es fundamental, pero sin ella no podemos calcular la distancia Unifrac. ![](https://i.imgur.com/xenuWs8.jpg) ::: Para pasar los datos de Metaphlan a phyloseq usamos la función metaphlantophyloseq. :arrow_down: Este script tienes que copiarlo y cargarlo en tu entorno de R, es decir te lo copias en un script de R y lo ejecutas sin más, luego más adelante lo usarás... ```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) } ``` ## Convertir en Objeto phyloseq * Metadata: Añadir metadata ojo con los Sample_names :eyes: Aqui un ejemplo con mis datos: ```r=1 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 ``` * Taxonomía: Leemos la tabla de abundancia "merged_abundance_table.txt". Quitamos la primera línea y arreglamos SampleID para que coincida. > Revisar que no haya errores con los números. Aqui un ejemplo: ```r=1 Abundance <- read_csv("Data/Metagenomica/merged_abundance_table.txt") Abundance <- column_to_rownames(Abundance, var="clade_name") Abundance$NCBI_tax_id <- NULL ``` Otro ejemplo: ```r=1 abundances.metahit<- read.table('metaphlan_phyloseq/merged_abundance_table.txt', skip = 1, sep = '\t', header = TRUE) rownames(abundances.metahit) <- abundances.metahit$clade_name abundances.metahit <- abundances.metahit %>% select(-c(clade_name)) abundances.metahit$NCBI_tax_id <- NULL ``` * Cargar script Function metaphlantophyloseq ```r=1 phyloseqin= metaphlanToPhyloseq(Abundance, metadat = sample) phyloseqin ``` * Añadir árbol filogenético para poder calcular distancia Unifrac ```r=1 #Arbol random_tree=rtree(ntaxa(phyloseqin), rooted = T, tip.label = taxa_names(phyloseqin)) #Objeto phyloseq con arbol physeq.tree <- merge_phyloseq(phyloseqin, random_tree) ``` ## Ver abundancias * Phylum más abundantes, puedes cambiarlo con los niveles taxonómicos que quieras, sirve igual para 16S o WGS. ```R library(tidyverse) #Porcentaje Phylum Phyl <- tax_glom(pseq.compositional, taxrank = "Phylum") tb_phyl <- psmelt(Phyl) tbphylum <- tb_phyl %>% group_by(Phylum) %>% summarise(M=median(Abundance)) %>% ungroup() %>% unique() %>% top_n(M,n=5) tbphylum <- tbphylum %>% inner_join(tb_phyl) Plotf <- ggplot(tbphylum, aes(SampleID, Abundance ,fill=Phylum)) + geom_col(position="fill") + scale_fill_manual(values=paleta5) Plotf <- Plotf + theme_minimal() + xlab("SampleID") + 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 ``` ## Análisis de diversidades: :::info El análisis de diversidad puedes hacerlo con `phyloseq` o con `microbiome` ::: > 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. ```r=1 GP1 <- transform_sample_counts(physeq.tree, function(x) 1E6 * x) ``` * Alfa diversidad ```r=1 plot_richness(GP1) # phyloseq rich <- richness(GP1) # microbiome plot_richness(GP1, 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() # phyloseq ``` * Transformamos datos composicionales con el paquete microbiome. ```r=1 library(microbiome) pseq.compositional <- transform(physeq.tree, "compositional") ``` * Beta diversidad ```r=1 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() ``` ## PICRUST con datos de 16S :pie: PICRUSt (pronounced “pie crust”) is a bioinformatics software package designed to predict metagenome functional content from marker gene (e.g., 16S rRNA) surveys and full genomes. Vamos a usar el plugin de QIIME2 para ejecutar PICRUSt [(q2-picrust2)](https://github.com/picrust/q2-picrust2) ### Instalación Lo saqué de [esta](https://library.qiime2.org/plugins/q2-picrust2/13/) página, ojo con las versiones de qiime2 y picrust, asegurate de usar las que ya tienes instaladas. ``` conda activate qiime2-amplicon-2024.5 conda install q2-picrust2=2024.5 \ -c conda-forge \ -c bioconda \ -c picrust ``` ### Ejecución Saqué este ejemplo del tutorial de esta [web](https://github.com/picrust/picrust2/wiki/q2-picrust2-Tutorial). Mejor echale un ojo y cambia por tus datos adecuadamente. ``` qiime picrust2 full-pipeline \ --i-table mammal_biom.qza \ --i-seq mammal_seqs.qza \ --output-dir q2-picrust2_output \ --p-placement-tool sepp \ --p-threads 1 \ --p-hsp-method pic \ --p-max-nsti 2 \ --verbose ``` **Con esto tendrías tes archivos:** `ec_metagenome.qza `- EC metagenome predictions (rows are EC numbers and columns are samples).Enzimas, aquí podriamos centrarnos en hidrolasas particulares que degraden las algas :turtle: `ko_metagenome.qza `- KO metagenome predictions (rows are KOs and columns are samples). Términos KO de vías metabolicas de [KEGG.](https://www.genome.jp/kegg/ko.html) `pathway_abundance.qza`- MetaCyc pathway abundance predictions (rows are pathways and columns are samples). Vías metabólicas de [metacyc. ](https://metacyc.org/) Quizás el archivo más interesante, ya que podemos ver vías metabolicas microbianas que posiblemente esten en las muestras tomando en cuenta las bacterias que hay presentes. Luego los archivos los pasamos a R y hacemos también un análisis exploratorio de los datos (abundancias, comparaciones a ver que tal). Y podemos hacer un análisis más especifico con [STAMP](https://beikolab.cs.dal.ca/software/STAMP). Y obtener este tipo de gráficos: ![image](https://hackmd.io/_uploads/HJgMd9fA0.png) ## LEfSe con datos de 16S LEfSe (Linear discriminant analysis Effect Size) determines the features (organisms, clades, operational taxonomic units, genes, or functions) most likely to explain differences between classes by coupling standard tests for statistical significance with additional tests encoding biological consistency and effect relevance. Se puede ejecutar desde los srevidores de [Galaxy](http://galaxy.biobakery.org/), pero los resultados son menos modificables que si los hacemos desde la línea de comandos. Te pongo el script para hacerlo desde CLI y si no podemos vamos a lo "fácil" con Galaxy :stuck_out_tongue_winking_eye: Este es mi codigo para pasar el objeto phyloseq con la estructura que necesita Lefse, hay que revisar que este todo correcto y a lo mejor tienes que hacer alguna modificación con tus datos: ```R=1 #Phyloseq to LEfSe stringsAsFactors = FALSE library(tidyverse) library(phyloseq) library(janitor) tax <- as.data.frame(t(physeq@otu_table@.Data)) tax <- rownames_to_column(tax, var="ASV") names <- as.data.frame(physeq@tax_table@.Data) names <- names %>% unite(Kingdom,Phylum,Class,Order,Family,Genus,col=Genus,sep="|") names <- rownames_to_column(names, var="ASV") tax <- names %>% select(ASV,Genus) %>% left_join(tax) tax$ASV <- NULL tax <- t(tax) tax <- tax %>% row_to_names(row_number=1) tax <- as.data.frame(tax) tax <- rownames_to_column(tax, var="SampleID") #Metadata metadata <- as.data.frame(physeq@sam_data) write_csv(metadata,"metadata.csv") #Select classes metadata <- read_csv("Juliana/metadata.csv") tax <- metadata %>% left_join(tax) lefse <- t(tax) lefse[is.na(lefse)] <- 0 lefse <-rownames_to_column(as.data.frame(lefse)) #Export write_tsv(lefse, "lefse.txt",col_names = F) ``` El código para ejecutar LEfSe en la línea de comandos: Primero instalar ```bash #Install with Conda conda install -c biobakery lefse ``` Ejecutar: ```bash=1 #Format data to LEfSe lefse_format_input.py lefse.txt lefse.in -c 2 -u 1 -s -1 -o 1000000 #-u SampleID -c Class -s Subclass (if there is not sublclass set to -1) # -o 1000000 scales the feature such that the sum (of the same taxonomic level) is 1M #Run LEfSe lefse_run.py lefse.in lefse.res -l 4 #-l Change LDA treshold (LDA > 2 by default) #Plot Results into bar chart lefse_plot_res.py lefse.res lefse.bar.png --format png --dpi 300 --title "" --class_legend_font_size 9 --width 10 #Plot cladogram lefse_plot_cladogram.py lefse.res lefse.cladogram.png --format png --dpi 300 --expand_void_lev 5 \ --class_legend_vis 0 --siblings_connector_width 0.5 --max_point_size 4 --alpha 0.05 \ --title "" --label_font_size 7 --labeled_stop_lev 7 --parents_connector_width 0.5 --abrv_stop_lev 6 ``` Cualquier cosa mejor ver el tutorial [aquí](https://huttenhower.sph.harvard.edu/lefse/).