###### tags: `16S rRNA` `WGS` `ML` `tutorials` # Para María [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/). * Análisis varios siguiendo el libro de [Orchestrating Microbiome Analysis.](https://microbiome.github.io/OMA/docs/devel/) ## Análisis 16S rRNA con qiime2 * Instalar qiime2 environment * Crear "manifest" en qiime2 * Visualizar calidad de las secuencias para elegir parametros "trim" y "trunc". * Ejecutar DADA2 para demux y quitar quimeras. * Construcción del árbol filogenético y asignación taxonómica * Calcular diversidades. ### Instalar qiime2 environment ```bash=1 # Descargar la última versión de qiime2 wget https://data.qiime2.org/distro/core/qiime2-2022.2-py38-linux-conda.yml # Instalar la última versión de qiime: conda env create -n <env_name> --file qiime2-2022.2-py38-linux-conda.yml # Optional clean-up rm qiime2-2022.2-py38-linux-conda.yml ``` :::success :bulb: qiime2 se actualiza a menudo, el archivo .yml cambiará en <qiime2-2022.2> a la versión más actual. Hoy 21/06/2023 <qiime2-2023.5>. ::: ### Crear "manifest" en qiime2 Para el paso posterior (calidad de las secuencias en qiime2), necesitamos crear primero un archivo "manifest.csv" que qiime2 utilizará para buscar los archivos con las secuencias. Para ello podemos crear un script que nos facilite el trabajo. ### Visualizar calidad de las secuencias Ya tenemos el fichero manifest.csv Tras esto, ya podemos importar los archivos en qiime2. :::success :bulb: Hasta ahora en todo momento estamos partiendo de raw data no procesado con anterioridad. Estos pasos son previos a eliminar los barcodes y primers. ::: A continuación, activamos el qiime environment creado e importamos los datos en qiime y, con *qiime demux summarize*, obtenemos un informe de la distribución de las calidades de las secuencias en cada posición: ```bash=1 conda activate <env_name> qiime tools import \ --type 'SampleData[PairedEndSequencesWithQuality]' \ --input-path manifest.csv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33 \ qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv ``` Tras ejecutar este código, obtenemos un fichero llamado "demux.qzv", el cual podemos visualizar. Para ello, podríamos usar el siguiente comando: ``` qiime tools view demux.qzv ``` Tras esto, he subido el fichero [aquí](https://view.qiime2.org/), obteniendo los siguientes resultados: ![](https://i.imgur.com/AFt838c.png) En la imagen se ve la calidad (eje Y) de las secuencias así como su tamaño en bases (eje X) para el forward y el reverse. Para entenderlo, ver ([este ejemplo](https://view.qiime2.org/visualization/?type=html&src=https%3A%2F%2Fdocs.qiime2.org%2F2022.2%2Fdata%2Ftutorials%2Fmoving-pictures%2Fdemux.qzv)). En él se puede ver sólo la calidad de las secuencias forward. Se ven secuencias muy pequeñas a la izq con buena calidad (**los primers**) y grandes de baja calidad a la dch (**secuencias alteradas o dañadas**). En comparación, nuestras secuencias están bastante limpias y con buena calidad en general. Es necesario determinar la calidad de las muestras para poder escoger en el siguiente paso los parámetros *p-trunc* (a partir de que número de bases la calidad no es buena) y *p-trim* (tamaño en bases de los primers). :::success :bulb: *p-trim* dependerá del tamaño de los primers: **Primers IlluminaMiseq V3-V4** –p-trim-left-f 19 –p-trim-left-r 22 **Primers IlluminaMiseq V4** –p-trim-left-f 19 –p-trim-left-r 19 ::: ### Ejecutar DADA2 para demux y quitar quimeras. Ejecutamos DADA2, eliminando los primers y las regiones de baja calidad de las secuencias que hemos observado. ```bash=1 qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-trim-left-f 19 \ # ponemos los valores de los primers establecidos --p-trim-left-r 22 \ --p-trunc-len-f 0 \ # Al estar las secuencias limpias dejamos los valores de trunc en 0, no dejarlo vacío --p-trunc-len-r 0 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qza ``` Aunque ya hayan sido eliminados los primers y las lecturas de baja calidad, hay que hacer el paso de DADA2 porque genera un fichero llamado *table.qza*, necesario para el siguiente paso en qiime2 (feature-table), y este paso a su vez también es necesario para la generación del árbol filogenético. Además, *table.qza* también es necesario para *phyloseq*. :::success :hourglass_flowing_sand: **Este paso lleva tiempo, dejarlo corriendo en sesión *tmux*** ::: Se generan 3 archivos: *rep-seqs.qza*, *denoising-stats.qza* y *table.qza*. Los usamos en el paso siguiente para generar las visualizaciones (archivos qzv). :eyes: La visualización de *table.qza* necesita que generemos un archivo *metadata.tsv* previamente. **Ese archivo tiene que tener al menos 2 columnas**: 1) '#SampleID' --> Se tiene que llamar así. Identifica las muestras. Si hay Fwd y Rev con que aparezca el ID solo una vez es suficiente. 2) Puede ser genotype, o sex, no hay un formato requerido. ```bash=1 qiime feature-table tabulate-seqs \ --i-data rep-seqs.qza \ --o-visualization rep-seqs.qzv qiime metadata tabulate \ --m-input-file denoising-stats.qza \ --o-visualization denoising-stats.qzv qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file metadata.tsv ``` **VISUALIZACIÓN DE LOS ARCHIVOS** (en este [link](https://view.qiime2.org/)) <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero rep-seqs.qzv</i></b></p> </font> ![](https://i.imgur.com/40kjoMo.png) Este fichero nos muestra una tabla con los IDs de las características o *features* (taxas) y sus secuencias asociadas ([fuente](https://docs.qiime2.org/2022.2/tutorials/moving-pictures/)). <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero denoising-stats.qzv</i></b></p> </font> ![](https://i.imgur.com/qaD9oD1.png) Este fichero nos enseña, para cada muestra, la variación en el número de secuencias tras los sucesivos pasos de filtrado llevados a cabo anteriormente ([fuente](https://forum.qiime2.org/t/interpreting-denoising-stats/6280/2)). Se incluye aquí el quitar las [secuencias quimera](https://learnmetabarcoding.github.io/LearnMetabarcoding/filtering/chimera_filtering.html). <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero table.qzv</i></b></p> </font> ![](https://i.imgur.com/OriBUYU.png) Este fichero nos da información acerca del número de muestras, el número de *taxa* (o *features*) existentes, el número de *taxa* identificados por muestra o el número de veces que una *taxa* se ha identificado. ### Construcción del árbol filogenético y asignación taxonómica Para la construcción del árbol filogenético, vamos usar la siguiente *pipeline*, que recurre al programa *mafft*: ```bash=1 qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences rep-seqs.qza \ --o-alignment aligned_filtered-rep-seqs.qza \ --o-masked-alignment masked-aligned-rep-seqs.qza \ --o-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza ``` :::success :clock1: **Para 20 muestras ~ 1h** ::: Los archivos generados en este paso serán necesarios para los comandos de *phyloseq*, concretamente el del árbol enraizado (*rooted-tree.qza*). ~~Nótese, además, que se ha tomado como *input* el fichero *filtered_rep-seqs.qza*, el cual ya ha sido filtrado respecto a las secuencias de mitocondrias y cloroplastos.~~ A continuación, vamos a llevar a cabo la **asignación taxonómica** con SILVA, lo cual es necesario hacer antes de los análisis de alfa y beta diversidad ya que se requiere para dichos pasos un fichero llamado *taxonomy.qza* que se genera con SILVA. Para ello, hay que descargar este clasificador de la página de qiime2 (bajar la versión más reciente): ```bash=1 wget https://data.qiime2.org/2022.2/common/silva-138-99-nb-classifier.qza ``` A continuación, ejecutamos qiime2 para la asignación taxonómica (con el fin de conocer la composición taxonómica de las muestras) usando este clasificador de tipo Naïve-Bayes ( lleva tiempo :clock1: :hourglass_flowing_sand:): ```bash=1 qiime feature-classifier classify-sklearn \ --i-classifier silva-138-99-nb-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv ``` Tras ejecutar estos comandos, obtenemos un fichero .qzv que puede ser visualizado [aquí](https://view.qiime2.org/), observándose lo siguiente: ![](https://i.imgur.com/EK4DSuF.png) Tal y como podemos apreciar, a cada taxón identificado se le asigna un valor de confianza. La mayoría de los taxones están expresados a nivel de género, si bien hay algunos en los que se especifica también la especie. Los taxones con mayor confianza están asociados a mitocondrias/cloroplastos, los cuales han aparecido aquí al poseer estos orgánulos también rRNA 16S. Tal y como discutimos en el GM, **vamos a eliminar estos taxones** que no son de relevancia aquí: ``` qiime taxa filter-table \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --p-exclude mitochondria,chloroplast \ --o-filtered-table filtered_table.qza qiime taxa filter-seqs \ --i-sequences rep-seqs.qza \ --i-taxonomy taxonomy.qza \ --p-exclude mitochondria,chloroplast \ --o-filtered-sequences filtered_rep-seqs.qza ``` Filtrada la **tabla** original, así como las **secuencias representativas**, podemos pedir a qiime que nos dibuje un barplot donde se representen las abundancias de cada ASV para cada caso (aquí he tenido que usar otro fichero llamado metadata.tsv (los datos separados por tabulaciones)): ``` qiime taxa barplot --i-table filtered_table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file datos_IAS.tsv \ --o-visualization taxa_barplot.qzv ``` Si lo visualizamos, obtenemos esto a nivel general, si bien podemos modificar los parámetros según nuestros intereses: <font color = 'steelblue'> <p style="text-align:center;"><b><i>Fichero taxa_barplot.qzv</i></b></p> </font> ![](https://i.imgur.com/ZzAbVX9.png) Tras obtener el fichero *taxonomy.qza*, así como *filtered_table.qza* y *filtered_rep-seqs.qza*, podemos realizar el análisis de alfa/beta-diversidad en *phyloseq*. ## Análisis en R :computer: :::info `Phyloseq` es un paquete de R que se usa mucho para análisis de microbioma. 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) ::: ### Análisis de Abundancias relativas A continuación, vamos a llevar a cabo el análisis de las abundancias relativas, así como de las alfa y beta diversidades por medio del paquete *phyloseq* de R. Para ello, ejecutamos R en la consola de Ubuntu e instalamos el paquete, además de otros paquetes necesarios para importar los datos: ```R=1 R #ejecutar R if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("phyloseq") #instalación de phyloseq if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("microbiome") #instalación de microbiome install.packages("remotes") remotes::install_github("jbisanz/qiime2R") #instalación de qiime2R install.packages("readr") #instalación de readr install.packages("tidyr") #instalación de tidyr install.packages("dplyr") #instalación de dplyr install.packages('knitr', dependencies = TRUE) #instalación de knitr conda install -c r r-nloptr # fuera de R install.packages("car") #instalación de car ``` Tras 15 minutos instalándose aproximadamente todos los paquetes, importamos las librerías: ``` library(phyloseq) library(tibble) library(qiime2R) library(tidyr) library(readr) library(dplyr) library(microbiome) library(knitr) library(car) library(stats) ``` Tras esto, pasamos los datos a la terminal de R para posteriormente poder hacer los análisis. Nótese el empleo de la **tabla ya filtrada** para generar la variable *ASV*: ```r=1 #Tabla de OTUs/ASV ASV<-read_qza("./filtered_table.qza") #Taxonomía taxonomyq <-read_qza("./taxonomy.qza") #Transformar tabla taxtableq <- taxonomyq$data %>% as_tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species")) taxtableq$Kingdom <- gsub("d__", "",taxtableq$Kingdom) taxtableq$Phylum <- gsub("p__", "",taxtableq$Phylum) taxtableq$Class <- gsub("c__", "",taxtableq$Class) taxtableq$Order <- gsub("o__", "",taxtableq$Order) taxtableq$Family <- gsub("f__", "",taxtableq$Family) taxtableq$Genus <- gsub("g__", "",taxtableq$Genus) taxtableq$Species <- gsub("s__", "",taxtableq$Species) #Clean taxa: Remove Kingdom Unassigned & Eukaryota taxtableq <- taxtableq[grep("Eukaryota|Unassigned", taxtableq$Kingdom, invert=T),] #Árbol filogenético que generamos anteriormente tree<-read_qza("./rooted-tree.qza") #Metadata metadata<-read_csv("./metadata.csv") #Or metadata<-read.delim("./metadata.tsv", sep="\t") #Construir Objeto phyloseq OTUs <- otu_table(ASV$data, taxa_are_rows = T) tree <- phy_tree(tree$data) TAXq <- tax_table(as.data.frame(taxtableq) %>% select_("-Confidence") %>% column_to_rownames("Feature.ID") %>% as.matrix()) #moving the taxonomy to the way phyloseq wants it sample_metadata <- sample_data(metadata %>% as.data.frame()) sample_names(sample_metadata) <- paste(metadata$ID) #Instead of ID put whatever is the name in metadata df ID column physeq<- merge_phyloseq(OTUs,tree,TAXq,sample_metadata) ``` Tras esto, he obtenido un objeto llamado *physeq* con el siguiente aspecto: ![](https://i.imgur.com/qbVexdP.png) ### Análisis de diversidades :::info El análisis de diversidad puedes hacerlo con `phyloseq` o con `microbiome` ::: * Alfa diversidad ```r=1 plot_richness(physeq) # phyloseq rich <- richness(physeq) # microbiome plot_richness(physeq, 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 ``` * Beta diversidad ```r=1 ordinated_taxa_unifrac = ordinate(physeq, method="MDS", distance="unifrac",weighted=TRUE) plot_ordination(physeq, ordinated_taxa_unifrac, color="Cluster", title = "Unifrac MDS Analysis") + theme_bw() ``` * Transformamos datos composicionales con el paquete microbiome. ```r=1 physeq.compositional <- transform(physeq, "compositional") ``` A partir de aquí ya puedes hacer el análisis que necesites dependiendo del objetivo del estudio, recomiendo [este](https://microbiome.github.io/OMA/docs/devel/) tutorial.