Para Adrian - HackMD

###### tags: `16S rRNA` `WGS` `ML` `tutorials` # Para Adrian [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 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 ellos podemos usar cualquiera de los siguientes scripts, según el lenguaje que prefiramos, creándolo o copiándolo en el directorio donde están los archivos. **ESCOGER UNO** ==**Perl**== Primero creamos el archivo .pl ```bash=1 touch script.pl #para crear el fichero vi script.pl #para abrir el fichero en la terminal y poder editarlo, puede ser vi, vim o nano ``` Una vez abierto en el editor copiar el siguiente código ```perl=1 #!/usr/bin/perl use Cwd 'abs_path'; open(OUT,">manifest.csv"); foreach $ar (@ARGV) { #Para sacar el $name del archivo: $name=$ar; @field = split (/_/, $ar); $name2 = @field[0]; #Para sacar el $path: $path =abs_path($ar); if ($name =~ /1.fq.gz/) { print OUT "$name2,$path,forward\n"; }else{ print OUT "$name2,$path,reverse\n"; } } close OUT; ``` Ejecutamos el script ```bash=1 perl script.pl *.fq.gz `` El archivo .csv generado estaría incompleto. Hay que añadirle la primera línea. Para ello: ```bash=1 vi manifest.csv # Para abrir el .csv sample-id,absolute-filepath,direction # Añadir como fila inicial en el archivo .csv # Salir y guardar, la fila se habrá añadido ``` ==**Python**== Primero creamos el archivo .py ```bash=1 touch 16S_manifest.py # Crear script vi 16S_manifest.py # abrirlo en el editor ``` Copiamos el siguiente código y salimos guardando. ```python=1 #!/usr/bin/python import os import csv cwd = os.getcwd() folders = [f for f in os.listdir() if not f.startswith('.')] folders = [f for f in folders if os.path.isdir(cwd+'/'+f)] with open('manifest.csv', 'w') as file: writer = csv.writer(file) writer.writerow(["sample-id", "absolute-filepath", "direction"]) for i in folders: os.chdir('./' + i) files = [f for f in os.listdir() if not f.startswith('.')] for j in files: if 'raw_1' in j: x = j.split('.') name = x[0] path = os.getcwd() + '/' + j direction = 'forward' writer.writerow([name, path, direction]) elif 'raw_2' in j: x = j.split('.') name = x[0] path = os.getcwd() + '/' + j direction = 'reverse' writer.writerow([name, path, direction]) else: continue os.chdir('../') ``` Ejecutamos el script. En este caso el archivo 'manifest.csv' que se genera ya está completo, no hay que añadir nada más. ```bash=1 python 16S_manifest.py ``` ### 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 ``` No es posible al estar en Gandalf.Copiar el archivo desde Gandalf hasta el PC local: ``` scp amartin@172.16.202.45:~/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/)) Fichero rep-seqs.qzv ![](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/)). Fichero denoising-stats.qzv ![](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). Fichero table.qzv ![](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: Fichero taxa_barplot.qzv ![](https://i.imgur.com/ZzAbVX9.png) Tras obtener el fichero *taxonomy.qza*, así como *filtered_table.qza* y *filtered_rep-seqs.qza*, podemos realizar el análisis de alfa/beta-diversidad en *phyloseq*. * **ANÁLISIS DE LAS ABUNDANCIAS RELATIVAS** A continuación, vamos a llevar a cabo el análisis de las abundancias relativas, así como de las alfa y beta diversidades por medio del paquete *phyloseq* de R. Para ello, ejecutamos R en la consola de Ubuntu e instalamos el paquete, además de otros paquetes necesarios para importar los datos: ```R=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*: ``` #Tabla de OTUs/ASV ASV<-read_qza("./filtered_table.qza") #Taxonomía taxonomyq <-read_qza("./taxonomy.qza") #Transformar tabla taxtableq <- taxonomyq$data %>% as_tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species")) taxtableq$Kingdom <- gsub("d__", "",taxtableq$Kingdom) taxtableq$Phylum <- gsub("p__", "",taxtableq$Phylum) taxtableq$Class <- gsub("c__", "",taxtableq$Class) taxtableq$Order <- gsub("o__", "",taxtableq$Order) taxtableq$Family <- gsub("f__", "",taxtableq$Family) taxtableq$Genus <- gsub("g__", "",taxtableq$Genus) taxtableq$Species <- gsub("s__", "",taxtableq$Species) #Árbol filogenético que generamos anteriormente tree<-read_qza("./rooted-tree.qza") #Metadata metadata<-read_csv("./metadata.csv") #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 metagenomas WGS. - Con notas de Blanca 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... 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): --- Comentarios de Blanca #### 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 ``` Fin comentarios de Blanca ------- ### Mis scripts y cosas #### MetaWRAP Sacar data de c/carpeta y juntar en una. Luego descomprimir Cambiar nombre/ quitar toda la :shit: del nombre > De D1807_ESFP200009043-1a_HCYJCDSXY_L2_1.fq a D1807_1.fq ```bash=v # movemos todos los fastq a carpeta reads: fastqfiles=$(find -mindepth 2 -type f | grep -v md5) for file in $fastqfiles; do mv $file reads; done ls reads | wc -l # todo ok, hay 250 (264 - 14 malos) # descomprimimos: cd reads gunzip *.gz ``` > - Para PRJEB5024: ```bash! # movemos todos los fastq a carpeta RAW_READS: fastqfiles=$(find -mindepth 2 -type f | grep -v md5 | grep '_') for file in $fastqfiles; do mv $file ../RAW_READS; done ls ../RAW_READS | wc -l # todo ok, hay 356: 178 runs descargados x 2 (forward + reverse) (178 porque partiamos de 186 bgi-MH - 8 con 'lane' de 06 y 12) # descomprimimos: cd ../RAW_READS gunzip *.gz ``` - Limpiar y filtrar reads del host... (Paciencia :hourglass_flowing_sand:) ```bash for F in RAW_READS/*_1.fq; do R=${F%_*}_2.fq BASE=${F##*/} SAMPLE=${BASE%_*} metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE & done ``` :::success ```bash! export PATH=$PATH:/home/lmarcos/metaWRAP/bin for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_03-11-2021-10-40.txt 2>&1 for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_19-11-2021-18-23.txt 2>&1 ``` ::: ------- Comentarios de Blanca #### Ejecución del comando Antes de correr todo esto hay que __añadir metawrap al PATH__. No me atreví a cambiar el path de forma definitiva por si acaso la liaba con algo, que alguna vez con mi ordenador me ha pasado. El comando que he usado es: :::success ```export PATH=$PATH:/home/lmarcos/metaWRAP/bin``` ::: > ***Nota:*** al correr `which config-metawrap` en este nuevo entorno y revisar el fichero, parece que sigue apuntando correctamente al hg18. Modifiqué un poco el comando porque los ficheros se han descomprimido como `.fastq`, no como `.fq`. Además, **antes de correrlo creé el directorio `READ_QC` porque me daba problemas si no lo hacía.** El comando que puse era: ```bash! for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE & done > meta20211102.log ``` > _3 nov 10:20_ Esto no ha funcionado. Pongo una prueba pequeña otra vez, no vaya a ser que haya cambiado algo del environment sin darme cuenta: > ```bash! > for F in RAW_READS_prueba/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC_prueba/$SAMPLE & done > read_qcERROR.txt 2>&1 > ``` >Esta prueba sí funciona. Mirando en el github de metaWRAP hubo una persona a la que le pasó algo parecido y le dijeron que [a veces hay problemas con el "&" para paralelizar](https://github.com/bxlab/metaWRAP/issues/344). :::success Probamos sin `&` (y usamos **`2>&1`** para **redirigir el standard error además del standard output**): ```bash! for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_03-11-2021-10-40.txt 2>&1 for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_19-11-2021-18-23.txt 2>&1 ``` >Este comando corrió sin problemas. Estuvo hasta las 21:40 h del 06/11 --> unas 83 horas, entre 75 muestras = como **1 h 6 min por muestra** > > ___Healthy MetaHit:___ corre hasta el 22/11/2021 a las 7:20 h, salen 175 carpetas. No aparecen mensajes de error. Habíamos descartado los runs ERR011189_1, ERR011088_2 y ERR011229_2, así que nos sale para las dos reverse que no existen (y por tanto no se hizo su control de calidad). ::: - Mover las secuencias limpias y sin human-host ```bash 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 ``` ```bash for f in READ_QC/concatenated-runs/; do mv $f CLEAN_READS/; done cd CLEAN_READS/concatenated-runs mv * .. cd .. rm -r concatenated-runs ``` :::danger **IMPORTANTE:** no nombrar los ficheros concatenados con una \_, porque el script de Perl que usamos para metaphlan se va a hacer un lío. ::: :::warning Quitamos la primera "\_" de los nombres: ```bash conc_files=$(ls CLEAN_READS/ | grep -E '_._..fastq') for file in $conc_files; do new_name=$(sed 's/\_//' <(echo $file)) mv CLEAN_READS/$file CLEAN_READS/$new_name done ``` Me he equivocado y he corrido un comando en el que la línea de `mv` era mv `CLEAN_READS/$file CLEAN_READS$new_name`. Movemos los ficheros que se han generado a CLEAN_READS y solucionamos los nombres: ```bash clean_files=$(ls CLEAN_READS/ | grep CLEAN') for file in $clean_files; do new_name=$(sed 's/CLEAN_READS//' <(echo $file)) mv CLEAN_READS/$file CLEAN_READS/$new_name done ``` ::: ### Metaphlan :custard: :::danger :warning:Verificar la última actualización de Metaphlan en la web de Biobakery, para intentar usar siempre la versión actualizada. Algunos comandos pueden cambiar dependiendo de la versión, así que esto es sólo una guia. ::: #### Instalación: Creamos entorno de conda: ```bash conda create --name metaphlan conda activate metaphlan conda install -c bioconda metaphlan ``` Ahora tenemos que instalar también la base de datos: > - "_If you have installed MetaPhlAn using Anaconda, it is advised to install the database in a folder outside the Conda environment. To do this, run:_ > `metaphlan --install --bowtie2db <database folder>` >- _If you install the database in a different location, remember to run MetaPhlAn using `--bowtie2db <database folder>`!_" >`metaphlan --install --bowtie2db ~/blanca/bowtiedb` #### Comandos: **Nota:** he modificado la línea: ```bash! system("metaphlan $r1,$r2 --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt") ``` Para incluir lo de `--bowtie2db`: ```bash! system("metaphlan $r1,$r2 --bowtie2db ~/blanca/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt") ``` De manera que al final el script `metaphlan-script.pl` tiene esta pinta: :::success ```perl! #!/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 ~/blanca/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt"); } ``` ::: Me he atrevido a correrlo con todas las muestras a la vez (este script tiene que estar en la carpeta `CLEAN_READS`): ```bash! ./metaphlan-script.pl *_1.fastq > metaphlanOUT_08-11-2021-13-09.txt 2>&1 ``` >**Nota 1:** antes de correrlo me he ido a la sesión de screen blanca (`screen -r blanca`), he desactivado el entorno metawrap-2, he activado el de metaphlan y ahí ya sí he ejecutado el comando >**Nota 2:** nada más empezar a correr, en `CLEAN_READS` había 850 GB > **Nota 3:** A las 15:25, lleva 7 archivos: va a unos 20 min por muestra --> 20 * 65 = 1300 min (21 h 40 min) :::success El script para las muestras de sanos cambia un poco, añadimos el argumento `--read_min_len 40` (ver diario 22/nov): ```perl! #!/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 ~/blanca/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt --read_min_len 40"); } ``` ```bash! #./metaphlan-script.pl *_1.fastq > metaphlanOUT_22-11-2021-13-44.txt 2>&1 ./metaphlan-script.pl *_1.fastq > metaphlanOUT_22-11-2021-16-23.txt 2>&1 ``` ::: >- Sesión de screen `blanca-metaphlan` >- Como daba fallos de 'ZeroDivisionError', he cambiado la longitud mínima para que entren todas las lecturas (hay ficheros donde la longitud de la lectura original es de 44 bp, y otros donde la longitud de original es 75 bp pero se ha reducido mucho al hacer el trimming). :::warning _Notas:_ - No hace falta hacer nada especial ni para el script ni para los argumentos (`$./script.pl *_1.fastq`) - Sobre el comando de metaphlan: - La entrada `metaphlan $r1,$r2` está en el formato necesario para manejar metagenomas paired-end - El argumento `--bowtie2out` permite guardar la salida de bowtie2, para que si tenemos que volver a correr este proceso sea más rápido al ahorrar este paso. - Para ello usaríamos `--input_type bowtie2out` (eso sí, tiene que estar generado con MetaPhlAn 3.0 en adelante, anteriores no son compatibles) ::: - Unir resultados: ```bash! python utils/merge_metaphlan_tables.py profiled_*.txt > output/merged_abundance_table.txt ``` > El comando que usé es ligeramente distinto: > ```bash! > merge_metaphlan_tables.py profiled_*.txt > ../OUTPUT/merged_abundance_table.txt > ``` ## 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: ```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) } ``` :::warning Te dejo mis scripts del último proyecto que he estado analizando por si te son de utilidad... Eso sí son un poco :spaghetti: ::: ## Convertir en Objeto phyloseq * Metadata: Añadir metadata ojo con los Sample_names :eyes: que si no coinciden peta... ```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 ``` >Los archivos de metadatos que tengo son: >- Para las muestras de MetaHit: > - `ena-metadata-selection.csv` tiene una selección de metadata del ENA, con una fila por muestra (los run_accession de las muestras concatenadas están actualizados) > - `metadata-costea.csv` tiene la información de la tabla EV1 de Costea, importante porque tenemos ahí BMI, edad, género y diabetes. > - Podemos hacer el join con `ena_sample_accession` o con `internal_sample_name` >- Para las muestras de Karlsson en principio lo tenemos todo en `curatedMetagenomicData` > Mi script para preparar metadatos: ```r=1 data.metahit.1 <- read.csv('./Datasets/ena-metadata-selection.csv') data.metahit.2 <- read.csv('./Datasets/metadata-costea.csv') data.metahit.2$sample_accession <- data.metahit.2$ena_sample_accession # une en una sola tabla: data.metahit <- data.metahit.1 %>% left_join(data.metahit.2, by = 'sample_accession') # elimina columnas no interesantes: data.metahit <- data.metahit %>% select(-c(internal_sample_name, ena_study_accession, ena_sample_accession, m_otu_diversity, inserts_post_qc, inserts_aligned_to_spec_i)) ``` :::warning **Blanca:** en este punto añadí un par de columnas al dataframe para que sirvieran de etiquetas una vez construyamos el clasificador (y también para poder comparar los grupos más fácilmente cuando hagamos análisis metagenómicos). Tenía dudas de cómo era mejor hacerlo, así que al final lo que decidí fue: - Añadir una columna `mu` con valores binarios, `1` indica **metabólicamente insano** y `0` indica sano. Con esta columna haríamos la clasificación en dos grupos. - En la columna `group` tenemos cuatro posibles valores: - `0` para no obesos metabólicamente sanos (**MHNO**) - `1` para obesos metabólicamente sanos (**MHO**) - `2` para no obesos metabólicamente insanos (**MUNO**) - `3` para obesos metabólicamente insanos (**MUO**) En este dataset de MetaHit, tendríamos solo 1's en la columna `mu`, mientras que en la columna `group` tendremos 2's y 3's. ::: ```r=1 # crea labels: data.metahit$mu <- rep(1, dim(data.metahit)[1]) data.metahit$group <- rep(5, dim(data.metahit)[1]) data.metahit$group[data.metahit$bmi < 30] <- 2 data.metahit$group[data.metahit$bmi >= 30] <- 3 # termina de arreglar el dataframe: data.metahit$group <- factor(data.metahit$group) data.metahit$mu <- factor(data.metahit$mu) row.names(data.metahit) <- data.metahit$run_accession # y aqui ya entra phyloseq: sample.metahit <- sample_data(data.metahit) sample_names(sample.metahit) <- gsub('E', 'profiled_E', data.metahit$run_accession) ``` * **Taxonomía:** Leemos la tabla de abundancia "merged_abundance_table.txt". Quitamos la primera línea y arreglamso SampleID para que coincida. > Revisar que no haya errores con los números. ```r=1 Abundance <- read_csv("Data/Metagenomica/merged_abundance_table.txt") Abundance <- column_to_rownames(Abundance, var="clade_name") Abundance$NCBI_tax_id <- NULL ``` > Blanca: ```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) ``` :::success Lo ideal sería en este punto unir con los resultados de `curatedmetagenomicData` ::: ## 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 ``` > Aquí me sale un aviso: > ``` > Warning message: > In estimate_richness(physeq, split = TRUE, measures = measures) : > The data you have provided does not have any singletons. This is highly suspicious. Results of richness estimates (for example) are probably unreliable, or wrong, if you have already trimmed low-abundance taxa from the data. > We recommended that you find the un-trimmed data and retry. > ``` > De momento no le doy mucha importancia por lo que he visto googleando un poco ([aquí](https://forum.qiime2.org/t/singletons-and-diversity-richness-indices/2971/5) y [aquí](https://github.com/benjjneb/dada2/issues/214)). * Transformamos datos composicionales con el paquete microbiome. ```r=1 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() ```