# Tutoriales análisis datos 16S rRNA - Anna :smile: ###### tags: `Estudiantes` `16S rRNA` `tutorials` * 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). * [Phyloseq](hhttps://joey711.github.io/phyloseq/) * [Microbiome package](https://microbiome.github.io/tutorials/). ## Renombrar archivos Con rename ``` $ sudo apt install rename ``` > Con -n vemos que va a hacer, si es lo que quieres ejecutalo de nuevo sin -n > ``` rename -n 's/L001_R2_001.fastq.gz/.fastq.gz/' *.fastq.gz ``` > Si quieres quitar algo de en medio.. ``` rename -n 's/CC1_//' *.fastq.gz ``` ### Data set CC-HFHF [Descargar](https://drive.google.com/file/d/1lenLHt3c91frTkzHbRoN3wUfU8TRq7zo/view?usp=sharing) - Crear archivo de metadata. > CCX CC: Control X:número de ratón. > HFX HF:Dieta grasa X:número de ratón. # Usar dataset con qiime2 * Importar datos a qiime2 con el "manifest" * Visualizar calidad de las secuencias para elegir parametros "trim" y "trunc". * Ejecutar DADA2 para demux y quitar quimeras. * Visualizar (uniendo con metadata en formato .tsv) ver sampling depth y escoger parametros. * Calcular diversidades. * Crear árbol filogenético. * Asignar taxonomía usando SILVA classifier. Instalar última version ``` wget https://data.qiime2.org/distro/core/qiime2-2021.4-py38-linux-conda.yml conda env create -n qiime2-2021.4 --file qiime2-2021.4-py38-linux-conda.yml #Optional clean-up rm qiime2-2021.4-py38-linux-conda.yml ``` * Crear manifest.csv con mi script :blush: ```perl=1 #!/usr/bin/perl use Cwd 'abs_path'; open(OUT,">manifiest.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; ``` Ejecutar ``` perl script_import_qiime.pl *.fq.gz ``` Poner en la primera fila del manifest :eyes: > sample-id,absolute-filepath,direction * Importar secuencias a Qiime2 ``` conda activate qiime2-2021.4 qiime tools import \ --type 'SampleData[PairedEndSequencesWithQuality]' \ --input-path manifiest.csv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33 \ qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv ``` Esto lleva su tiempito :zzz: Para escoger los parámetros de *p-trunc* ver la calidad de las secuencias. Para los parámetros de *p-trim* ver tamaño de los primers y quitarlos. > 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 ``` qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-n-threads 6 \ --p-trim-left-f 19 \ --p-trim-left-r 22 \ --p-trunc-len-f 220 \ --p-trunc-len-r 220 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qza ``` * Preparar archivo de metadata. Formato .tsv indicando grupos y número de muestra. > Al menos tiene que haber dos columnas :eyes: y la primera línea debe empezar con: #SampleID ``` qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file metadata-PYM.tsv 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 ``` * Construir árbol filogenético :deciduous_tree: ``` qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences rep-seqs.qza \ --o-alignment aligned-rep-seqs.qza \ --o-masked-alignment masked-aligned-rep-seqs.qza \ --o-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza ``` * Diversidades con qiime (Yo prefiero después calcularlo con el phyloseq) El sampling depth lo escogemos en base a las visualizaciones previas. > En este caso he decidio que el sampling depth sea 50344 de esta forma matengo todas las muestras en el análisis y un 76% de las features observadas. ``` qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 50344 \ --m-metadata-file metadata-PYM.tsv \ --output-dir core-metrics-results ``` * **Asignación taxonómica** > :eyes: Para que funcione el clasificador debe ser la misma version del scikit-learn instalada con qiime. La última actualización de ambos coincide (Qiime y scikit-learn). * Descargar SILVA classifier (Actualizado) >GreenGenes está desactulizado desde 2013, por lo que usamos SILVA. ``` wget https://data.qiime2.org/2021.4/common/silva-138-99-nb-classifier.qza ``` * Asignación taxonómica :clock1: ``` 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 qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file metadata-PYM.tsv \ --o-visualization taxa-bar-plots.qzv ``` **Colapsar data a nivel de Género** >Algunas secuencias se asignan a nivel de especie. Pero sabemos que al usar 16S rRNA no somos capaces de tener tanta discriminación como para asignar especies, por lo que es preferible usar la opción de taxa-collapse. Esta función va a unir todos los "features" que tengan la misma asignación taxonómoca al nivel que digamos (En este caso 6 - género-). * Función taxa-collapse a nivel de género ``` qiime taxa collapse \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --p-level 6 \ --o-collapsed-table collapsed.qza ``` # Análisis en R * Librerías ``` library(tidyverse) library(qiime2R) library(phyloseq) library(microbiome) library(viridis) ``` * Pasar datos de qiime2 a phyloseq ``` #Tabla de OTUs/ASV ASV<-read_qza("Data/Qiime/table.qza") #Taxonomia taxonomyq <-read_qza("Data/Qiime/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) #Arbol tree<-read_qza("Data/Qiime/rooted-tree.qza") #Metadata metadata<-read_csv("Data/metadata_informe.csv") #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$SampleID) physeq<- merge_phyloseq(OTUs, tree, TAXq,sample_metadata) ``` ### Análisis de resultados * Diversidad ``` #Alfa diversidad rich <- alpha(physeq, index = "all") plot_richness(physeq, color = "Cluster", x = "Cluster", measures = c("Chao1", "Simpson", "Shannon")) + geom_boxplot(aes(fill = Cluster), alpha=.7) + theme_bw() pseq.rarified <- rarefy_even_depth(physeq) rich_rare <- alpha(pseq.rarified, index = "all") plot_richness(pseq.rarified,, color = "Cluster", x = "Cluster", measures = c("Chao1", "Simpson", "Shannon")) + geom_boxplot(aes(fill = Cluster), alpha=.7) + theme_bw() ``` * Estadísticos para ver diferencias en sólo los indices de Chao1 Shanon y Simpson ``` alpha <- estimate_richness(pseq.rarified, measures=c("Chao1", "Shannon", "Simpson")) alpha <- alpha %>% as.data.frame() %>% rownames_to_column("SampleID") %>% left_join(metadata) #Diferencias en varianzas library(car) leveneTest(y = alpha$Chao1, group = alpha$Grupo, center = "median") #Todo junto Chao1 <- kruskal.test(Chao1 ~ Grupo, data = alpha) Shannon <- kruskal.test(Shannon ~ Grupo, data = alpha) Simpson <- kruskal.test(Simpson ~ Grupo, data = alpha) #Wilcoxon en los que he encontrado dif. dif_Chao <- pairwise.wilcox.test(x = alpha$Chao1, g = alpha$Grupo, p.adjust.method = "bonferroni" ) ``` * Transformar datos ``` #Compositional pseq.compositional <- transform(physeq, "compositional") ``` * Beta-diversidad ``` #Unifrac ordinated_taxa_unifrac = ordinate(pseq.compositional, method="MDS", distance="unifrac",weighted=TRUE) plot_ordination(pseq.compositional, ordinated_taxa_unifrac, color="BMI", title = "Unifrac MDS Analysis") + theme_bw() ``` * Phylum más abundantes ``` #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 ``` ## LEfSe * Filtrar taxa con una abundancia <0.001. :eyes: Esto lo ajustamos para cada estudio dependiendo del número de ASV/OTUs que perdamos con el filtrado y somos más o menos restrictivos al respecto. ``` pseq.fil <- filter_taxa(pseq.compositional, function(x) mean(x) > 0.001, TRUE) ``` Transformar objeto Phyloseq para usar en [LEfSe](https://huttenhower.sph.harvard.edu/lefse/). El input para LEfSE tiene que ser >Sample_ID >Grupo >ASV1 >ASV2 conda install -n lefse -c biobakery lefse Exportar como .tsv ``` #Export write_tsv(lefse, "lefse.txt",col_names = F) ``` Aquí el [código](https://github.com/laurichi13/Phyloseq-to-LEfSe) que suelo usar. ## Diferencias entre bichos específicos Para ver diferencias de un bicho entre dos grupos (Esto sería válido sólo para la misma visita). Usar kruskal-wallis y si salen diferencias ver el Wilcoxon test. ``` #Lactococcus kruskal.test(Lactococcus ~ Grupo, data = dif_gen) pairwise.wilcox.test(x = dif_gen$Lactococcus, g = dif_gen$Grupo, p.adjust.method = "holm" ) Lacto <- ggplot(dif_gen,aes(Grupo, Lactococcus)) Lacto <- Lacto + geom_boxplot(aes(fill=Grupo)) + theme_light() + scale_fill_brewer(palette="Set1") Lacto <- Lacto + theme(axis.text.x = element_text(color = "black", size = 14, angle = 0, hjust = .5, vjust = .5, face = "plain")) + theme(axis.text.y = element_text(color = "black", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain")) + ggtitle ("Lactococcus sp. relative abundance") + theme(plot.title = element_text(family="Arial", size=rel(2), vjust=2,face="plain", color="black", lineheight=1.5)) Lacto <- Lacto + labs(x = "Grupo",y = "Relative abundance") + theme(axis.title = element_text(size=rel(1.5))) Lacto ``` ## Para muestras pareadas Ver este [tutorial](https://bookdown.org/dietrichson/metodos-cuantitativos/prueba-t-para-muestras-pareadas.html). Esto sería para ver como cambia un bicho entre las diferentes visitas. ## Otros análisis ### MFA análisis para escoger variables más importantes [Aqui ](http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/116-mfa-multiple-factor-analysis-in-r-essentials/)se explica en que consiste y como se hace. Un script mio. ``` library(tidyverse) library(FactoMineR) library(factoextra) DATA_plus <- read_csv("Data/Objetos/data_plus.csv") #Multiple Factor Analisys DATA_plus <- column_to_rownames(DATA_plus, var="Sample_ID") colnames(DATA_plus) #"s" para estandarizar variables continuas; "n" para variables categoricas res.mfa <- MFA(DATA_plus, group = c(10,11,1,3,3,1,1,3,1,1,10,1,15,1,1,1,1,5,7,8,1), type = c("s","s","s","s","s","s","s","s","n", "s", "n", "s","n","s","s","s","s","s","s","s","s"), name.group = c("hema","bioq","edad","antropo","grasa_musculo","grasa.visceral","metabolismo","presion.fc","CDAT","GIP.heces","sintomas1", "sumatoriaS1","sintomas2","sumatoriaS2a", "sumatoriaS2b", "CeD", "GSRS", "inflamacion","citoquinas","mucosa","IFAP"), num.group.sup = NULL, graph = FALSE) print(res.mfa) ``` ## Clustering > :eyes: Hay muchas técnicas de clustering! - Priemro ver [esto](http://www.sthda.com/english/wiki/wiki.php?id_contents=7932), para ver que técnica sería la más adecuada con tus datos - Aquí se explica como se hace y en que consiste el [MClust](https://cran.r-project.org/web/packages/mclust/vignettes/mclust.html), me gusta este método. - [Hierarchical clustering](https://www.datanovia.com/en/courses/hierarchical-clustering-in-r-the-essentials/). - [ k -means](https://www.datanovia.com/en/lessons/k-means-clustering-in-r-algorith-and-practical-examples/). # Análisis futuros - Clasificación de pacientes en base a marcadores clínicos e inlflamatorios. Hacer clustering y ver que pasa... sin tomar en cuenta respondedores vs. no respondedores. - k-means de abundancia y/o unifrac - Correlación entre muestras todo vs. todo. - Análisis tipo este artículo de [Eran Segal ](https://www.cell.com/cell-metabolism/fulltext/S1550-4131(17)30288-7?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS1550413117302887%3Fshowall%3Dtrue) y efecto del pan sobre marcadores metabólicos. i.e: Evaluar variabilidad individual.