--- title: "Chicken intestinal microbiota" author: "Daniel Guariz Pinheiro" format: html editor: visual bibliography: references.bib --- ```{r, label="loading-R-libraries"} #| echo: false #| eval: true # Carregando bibliotecas: suppressMessages( library("tidyverse") ) # Ciência de Dados em R suppressMessages( library("dada2")) # Processamento de ASVs suppressMessages( library("phyloseq") ) # Análises de microbioma suppressMessages( library('vegan') ) # Análises ecológicas suppressMessages( library('openxlsx') ) # Manipulação de arquivos do Excel suppressMessages( library("ShortRead") ) # Manipulação de arquivos de sequências suppressMessages( library('DECIPHER') ) suppressMessages( library('phangorn') ) suppressMessages( library('metagMisc') ) suppressMessages( library('microeco') ) suppressMessages( library('file2meco') ) suppressMessages( library('cowplot') ) suppressMessages( library('hilldiv') ) suppressMessages( library('phytools') ) suppressMessages( library('picante') ) suppressMessages( library('microbiome') ) suppressMessages( library('ggpubr') ) suppressMessages( library('randomcoloR') ) suppressMessages( library('UpSetR') ) suppressMessages( library('MiscMetabar') ) suppressMessages( library('ComplexUpset') ) suppressMessages( library('MicrobiotaProcess') ) suppressMessages( library('ComplexHeatmap') ) suppressMessages( library('RColorBrewer') ) suppressMessages( library('ggthemes') ) suppressMessages( library('GUniFrac') ) suppressMessages( library('DESeq2') ) suppressMessages( library('multcompView') ) suppressMessages( library('colorRamp2') ) suppressMessages( library('Rmisc') ) suppressMessages( library('mgcv') ) suppressMessages( library('vsn') ) suppressMessages( library('data.table' ) ) suppressMessages( library('NMF') ) ``` ```{r, label="loading-extra-functions"} #| echo: false #| eval: true # Função para contar ASVs por amostra getN <- function(x) sum(getUniques(x)) # Função: Resumo das taxonomias summary.tax <- function(ordered_names, ...){ tmp_names <- ordered_names tmp_list <- list(...) tmp_summary1 <- data.frame(Database = NA, Kingdom = NA, Phylum = NA, Class = NA, Order = NA, Family = NA, Genus = NA, Species = NA) for (i in seq(length(tmp_list))) { tmp_db <- tmp_list[[i]] %>% select(Kingdom:Species) tmp_name <- tmp_names[[i]] tmp_ksum <- sum(table(tmp_db["Kingdom"])) tmp_psum <- sum(table(tmp_db["Phylum"])) tmp_csum <- sum(table(tmp_db["Class"])) tmp_osum <- sum(table(tmp_db["Order"])) tmp_fsum <- sum(table(tmp_db["Family"])) tmp_gsum <- sum(table(tmp_db["Genus"])) tmp_ssum <- sum(table(tmp_db["Species"])) tmp_summary1 <- tmp_summary1 %>% add_row(Database = tmp_name, Kingdom = tmp_ksum, Phylum = tmp_psum, Class = tmp_csum, Order = tmp_osum, Family = tmp_fsum, Genus = tmp_gsum, Species = tmp_ssum) } tmp_summary1 <- tmp_summary1[-1,] tmp_summary2 <- data.frame(Database = NA, UnqKingdom = NA, UnqPhylum = NA, UnqClass = NA, UnqOrder = NA, UnqFamily = NA, UnqGenus = NA, UnqSpecies = NA) for (i in seq(length(tmp_list))) { tmp_db <- tmp_list[[i]] %>% select(Kingdom:Species) tmp_name <- tmp_names[[i]] tmp_ksum_u <- sum(table(unique(tmp_db["Kingdom"]))) tmp_psum_u <- sum(table(unique(tmp_db["Phylum"]))) tmp_csum_u <- sum(table(unique(tmp_db["Class"]))) tmp_osum_u <- sum(table(unique(tmp_db["Order"]))) tmp_fsum_u <- sum(table(unique(tmp_db["Family"]))) tmp_gsum_u <- sum(table(unique(tmp_db["Genus"]))) tmp_ssum_u <- sum(table(unique(tmp_db["Species"]))) tmp_summary2 <- tmp_summary2 %>% add_row(Database = tmp_name, UnqKingdom = tmp_ksum_u, UnqPhylum = tmp_psum_u, UnqClass = tmp_csum_u, UnqOrder = tmp_osum_u, UnqFamily = tmp_fsum_u, UnqGenus = tmp_gsum_u, UnqSpecies = tmp_ssum_u) } tmp_summary2 <- tmp_summary2[-1,] return(list(tmp_summary1, tmp_summary2)) } # Função: Parseamento das taxonomias fix_tax <- function(tax) { tmp <- as.data.frame(tax) %>% dplyr::mutate(across(Kingdom:Species, ~ str_replace_na(., "Unclassified"))) tmp[2] <- str_remove_all(tmp[[2]],"k__") tmp[3] <- str_remove_all(tmp[[3]],"p__") tmp[4] <- str_remove_all(tmp[[4]],"c__") tmp[5] <- str_remove_all(tmp[[5]],"o__") tmp[6] <- str_remove_all(tmp[[6]],"f__") tmp[7] <- str_remove_all(tmp[[7]],"g__") tmp[8] <- str_remove_all(tmp[[8]],"s__") tmp[2] <- paste0("k__",tmp[[2]]) tmp[3] <- paste0("p__",tmp[[3]]) tmp[4] <- paste0("c__",tmp[[4]]) tmp[5] <- paste0("o__",tmp[[5]]) tmp[6] <- paste0("f__",tmp[[6]]) tmp[7] <- paste0("g__",tmp[[7]]) tmp[8] <- paste0("s__",tmp[[8]]) names(tmp) <- c("ID", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species") tmp[grep("_Unclassified$", tmp$Species), 8] <- "s__Unclassified" tmp$Species <- sub("s__\\(.*", NA, tmp[,8]) tmp$Species <- gsub("\\(.*", "", tmp[,8]) return(tmp) } loessErrfun_mod4 <- function(trans) { qq <- as.numeric(colnames(trans)) est <- matrix(0, nrow=0, ncol=length(qq)) for(nti in c("A","C","G","T")) { for(ntj in c("A","C","G","T")) { if(nti != ntj) { errs <- trans[paste0(nti,"2",ntj),] tot <- colSums(trans[paste0(nti,"2",c("A","C","G","T")),]) rlogp <- log10((errs+1)/tot) # 1 psuedocount for each err, but if tot=0 will give NA rlogp[is.infinite(rlogp)] <- NA df <- data.frame(q=qq, errs=errs, tot=tot, rlogp=rlogp) # original # ###! mod.lo <- loess(rlogp ~ q, df, weights=errs) ###! # mod.lo <- loess(rlogp ~ q, df, weights=tot) ###! # # mod.lo <- loess(rlogp ~ q, df) # jonalim's solution # https://github.com/benjjneb/dada2/issues/938 mod.lo <- loess(rlogp ~ q, df, weights = log10(tot),degree = 1, span = 0.95) pred <- predict(mod.lo, qq) maxrli <- max(which(!is.na(pred))) minrli <- min(which(!is.na(pred))) pred[seq_along(pred)>maxrli] <- pred[[maxrli]] pred[seq_along(pred)MAX_ERROR_RATE] <- MAX_ERROR_RATE est[est% data.frame() %>% mutate_all(funs(case_when(. < X40 ~ X40, . >= X40 ~ .))) %>% as.matrix() rownames(est) <- rownames(estorig) colnames(est) <- colnames(estorig) # Expand the err matrix with the self-transition probs err <- rbind(1-colSums(est[1:3,]), est[1:3,], est[4,], 1-colSums(est[4:6,]), est[5:6,], est[7:8,], 1-colSums(est[7:9,]), est[9,], est[10:12,], 1-colSums(est[10:12,])) rownames(err) <- paste0(rep(c("A","C","G","T"), each=4), "2", c("A","C","G","T")) colnames(err) <- colnames(trans) # Return return(err) } ggrarecurve.get_rarecurve <- function(obj.ps, step=100, file=NULL){ set.seed(1024) if (!is.null(file)) { svg(filename=file) } else { pdf(NULL) } rcx <- vegan::rarecurve(t(as.matrix(as.data.frame(otu_table(obj.ps)))), step = step, label = T) graphics.off() names(rcx) <- colnames(otu_table(obj.ps)) samp_data <- sample_data(obj.ps) samp_data_names <- colnames(samp_data) rareres <- data.frame(matrix(nrow=0, ncol=length(samp_data_names)+3)) colnames(rareres) <- c(samp_data_names,'readsNums','value','sample') cnt<-0 for (n in names(rcx)) { for (ni in seq(1,length(names(rcx[[n]]))) ) { cnt<-cnt+1 rareres[cnt,'sample'] <- n for (cn in samp_data_names) { rareres[cnt,cn] <- samp_data[n,cn] } rareres[cnt,'readsNums'] <- attributes(rcx[[n]])$Subsample[ni] rareres[cnt,'value'] <- as.numeric(rcx[[n]][ni]) rareres[cnt,'Alpha'] <- 'Observe' } } return(rareres) } ggrarecurve.rarecurve <- function(obj, #sampleda, indexNames="Observe", linesize=0.5, facetnrow=1, shadow=TRUE, shadow.alpha=0.3, factorNames, factorColor=NULL, factorLine=NULL, se=FALSE, ...){ method="lm" formula=y ~log(x) if (is.null(factorColor)) { factorColor <- factorNames[1] } if (is.null(factorLine)) { factorLine <- factorNames[1] } mapping <- aes_string(x="readsNums", y="value", color="sample") if (!is.null(indexNames)){ obj <- obj %>% dplyr::filter(.data$Alpha %in% indexNames) } if (!missing(factorNames)){ if (shadow){ obj.sum <- summarySE(obj, measurevar="value", groupvars=c(factorNames, "readsNums", "Alpha"), .drop=TRUE, na.rm=TRUE) #obj.sum <- subset(obj.sum, ! is.na(se) ) obj.sum[which(obj.sum$se>obj.sum$value),'se'] <- (obj.sum[which(obj.sum$se>obj.sum$value),'value']/1.96) obj.sum$up <- obj.sum$value - (1.96*obj.sum$se) obj.sum$down <- obj.sum$value + (1.96*obj.sum$se) obj.sum[ which(obj.sum$up <= 0), 'up' ] <- 1 obj.ci <- data.frame(matrix(nrow=0, ncol=length( c(factorNames, 'readsNums', 'down', 'up') ))) colnames(obj.ci) <- c(factorNames, 'readsNums', 'down', 'up') for (f in unique(obj.sum[[ factorNames[1] ]]) ) { obj.tmp <- obj.sum %>% dplyr::filter(.data[[ factorNames[1] ]] == f) obj.tmp <- subset(obj.tmp, ! is.na(se) & N>=1) # fit formula m.down <- lm(down ~ log(readsNums), data = obj.tmp) m.up <- lm(up ~ log(readsNums), data = obj.tmp) # define finer grid of predictor values xnew <- seq(1,max(obj.tmp$readsNums), by = 10) # apply predict() function to the fitted LM model # using the finer grid of x values p.down <- predict(m.down, newdata = data.frame(readsNums = xnew), se = FALSE) p.up <- predict(m.up, newdata = data.frame(readsNums = xnew), se = FALSE) g <- data.frame(readsNums = xnew, down = p.down, up = p.up) for (cl in factorNames ) { g[[cl]] <- unique(obj.tmp[[cl]]) } obj.ci <- rbind(obj.ci,g) } obj.ci$Alpha <- 'Observe' mapping <- modifyList(mapping, aes_string(group=factorNames[1], color=factorColor, linetype=factorLine, fill=factorColor)) }else{ obj.sum <- obj mapping <- modifyList(mapping, aes_string(group="sample", color=factorColor, fill=factorColor, linetype=factorLine)) } } p <- ggplot(data=obj.sum, mapping=mapping) if (!missing(factorNames) && shadow){ p <- p + geom_ribbon(data=obj.ci, mapping=aes_string(fill=factorColor, ymin='up', ymax='down', x='readsNums', color=factorColor), inherit.aes=FALSE, alpha=shadow.alpha, color=NA, show.legend=FALSE) p <- p + geom_line(data=obj.ci, mapping=aes_string(y='up', x='readsNums', color=factorColor, linetype=factorLine), inherit.aes=FALSE, show.legend=FALSE,linewidth=0.01) p <- p + geom_line(data=obj.ci, mapping=aes_string(y='down', x='readsNums', color=factorColor, linetype=factorLine), inherit.aes=FALSE, show.legend=FALSE,linewidth=0.01) } p <- p + geom_smooth(se=se, method = method, linewidth=linesize,formula = formula, ...)+ scale_y_continuous(limits=c(0,NA), oob=scales::squish) + facet_wrap(~ Alpha, scales="free", nrow=facetnrow) + ylab("alpha metric")+xlab("number of reads") return(p) } setSharedPerc <- function(dt.perc.obj, phyloseq.obj, down.l=1, up.l=dim(sample_data(phyloseq.obj))[1] ) { library('phyloseq') library('tidyverse') library('data.table') proc.samps <- c() for (l.sample_i in down.l:up.l) { l.sample <- rownames(dt.perc.obj)[l.sample_i] for (c.sample_i in 1:dim(dt.perc.obj)[2]) { c.sample <- colnames(dt.perc.obj)[c.sample_i] #cat(paste0("CALL setSharedPerc(): ",down.l," and ",up.l," ",l.sample," ",c.sample,"=",dt.perc.obj[l.sample,c.sample],"\n")) if ( is.na(dt.perc.obj[l.sample_i, ..c.sample_i]) ) { #cat(paste0("Analyzing ",l.sample," vs ",c.sample,"\n")) n.l <- 0 n.c <- 0 n.l <- as.numeric(table(as.data.frame(otu_table(phyloseq.obj))[[l.sample]]>0)['TRUE']) n.c <- as.numeric(table(as.data.frame(otu_table(phyloseq.obj))[[c.sample]]>0)['TRUE']) n.l_c <- 0 if (l.sample != c.sample) { # newDF <- subset(as(sample_data(phyloseq.obj), "data.frame"), X.NAME %in% c(l.sample,c.sample)) # phylo.ss <- phyloseq.obj # sample_data(phylo.ss) <- sample_data(newDF) #phylo.ss <- subset_samples(phyloseq.obj, X.NAME %in% c(l.sample,c.sample) ) n.l_c <- as.numeric(table(apply( as.data.frame(otu_table(phyloseq.obj))[,c(l.sample,c.sample)], 1, function(x) sum(x >= 1) == (2)))['TRUE']) %>% replace_na(0) } else { n.l_c <- n.l } #dt.perc.obj[l.sample,c.sample] <- (n.l_c/n.l)*100 #dt.perc.obj[c.sample,l.sample] <- (n.l_c/n.c)*100 data.table::set( x = dt.perc.obj, i = l.sample_i, j = c.sample_i, value = (n.l_c/n.l)*100) data.table::set( x = dt.perc.obj, i = c.sample_i, j = l.sample_i, value = (n.l_c/n.c)*100) proc.samps <- unique(c(proc.samps,l.sample)) } } } return(proc.samps) } ``` ## [Amostras e objetivos]{style="color: #181C14"} Este estudo tem como objetivo compreender o microbioma de segmentos do intestino de frangos que ainda não eclodiram do ovo tratados com probióticos. A análise busca revelar as diferenças e semelhanças nas comunidades microbianas associadas a cada condição. ```{r, label="setup-R", include=FALSE} #| echo: false #| warnings: false #| eval: true #| include: false # Primeiros passos base.dir <- "/data/MicrobiomaFrangos" Sys.setenv(BASE_DIR=base.dir) primers.file <- paste0(base.dir,"/refs/primers.fa") primers <- readFasta(primers.file) Sys.setenv(PRIMERS_FILE=primers.file) ``` ## Pré-processamento ### Nomenclatura As amostras referentes às bibliotecas de *amplicons* 16S sequenciadas foram nomeadas de acordo com o padrão: | CODE | TREATMENT | |------|------------| | C | Control | | V | Vehicle | | P | Probiotic | | A | Antibiotic | | CODE | SEGMENT | |------|----------| | D | Duodenum | | J | Jejunum | | I | Ileum | | C | Cecum | | CODE | SEX | |------|--------| | M | Male | | F | Female | \\\\ ## Processamento dos dados ### Análise de qualidade das amostras pré-processamento A avaliação de qualidade foi realizada com FastQC v0.11.9 utilizando os parâmetros pré-definidos (*default*) para as leituras R1 e R2 de cada amostra, e, por fim, os relatórios individuais foram reunidos em um relatório final com MultiQC v1.13.dev0. Programas: FastQC: *A high throughput sequence QC analysis tool*. - fastqc v0.11.9 MultiQC: *Aggregate results from bioinformatics analyses across many samples into a single report*. - multiqc v1.13.dev0 ```{bash, label='fastqc-before', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/input" outdir="${BASE_DIR}/output/fastqc/pre" mkdir -p ${outdir} rm -f ${outdir}/R1.txt rm -f ${outdir}/R2.txt for fq in $(ls ${indir}/*_R1*.fastq) ; do bn=$(basename ${fq} _R1.fastq) echo "Running FastQC for ${bn}_R1 and ${bn}_R2 ..." fastqc \ ${indir}/${bn}_R1.fastq \ ${indir}/${bn}_R2.fastq \ --threads 10 \ --outdir ${outdir} > /dev/null 2> /dev/null echo "${outdir}/${bn}_R1_fastqc.zip" >> ${outdir}/R1.txt echo "${outdir}/${bn}_R2_fastqc.zip" >> ${outdir}/R2.txt done multiqc --title "Pre-processing reads (R1)" --file-list ${outdir}/R1.txt -o ${outdir} multiqc --title "Pre-processing reads (R2)" --file-list ${outdir}/R2.txt -o ${outdir} ``` #### Resultados da Avaliação de qualidade (prévia) [MultiQC Report - R1](./fastqc/pre/Pre-processing-reads-R1_multiqc_report.html) [MultiQC Report - R2](./fastqc/pre/Pre-processing-reads-R2_multiqc_report.html) ### Verificação dos iniciadores (*primers*) O módulo *search_oligodb* do programa usearch v11.0.667 foi utilizado para a busca pelas sequências dos pares de oligos iniciadores (*primers*) utilizados na amplificação dos fragmentos. Essa busca ocorreu considerando as primeiras 500 sequências, no máximo 2 *mismatches* e buscando em ambas as orientações. Programa: USEARCH: *search and clustering algorithms that are often orders of magnitude faster than BLAST*. - USEARCH v11.0.667, módulo *-search_oligodb*. ```{bash, label='usearch-search_oligodb', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/input" outdir="${BASE_DIR}/output/usearch/pre/findadapt" primers="${BASE_DIR}/refs/primers.fa" mkdir -p ${outdir} for fq in $(ls ${indir}/*.fastq) ; do bn=$(basename ${fq} .fastq) echo "Finding primers into file basename ${bn}" head -n 2000 ${fq} > ${outdir}/tmp.fq head -n 2000 ${fq} > ${outdir}/tmp.fq usearch11 \ -search_oligodb ${outdir}/tmp.fq \ -db ${primers} \ -maxdiffs 2 \ -threads 10 \ -strand both \ -userout ${outdir}/${bn}_findAdapts.txt \ -userfields query+target+qstrand+diffs+tlo+thi+qlor+qhir+trowdots > /dev/null 2> /dev/null rm -f ${outdir}/tmp.fq done ``` Os *primers* utilizados na verificação foram: ```{r, label="print-primers"} #| echo: false #| eval: true print(sread(primers)) ``` ### Resultado da verificação dos iniciadores (*primers*) [Arquivos](./usearch/pre/findadapt/) ### Estimativas de tamanhos e qualidades das leituras (*reads*) O módulo *fastq_eestats2* do programa usearch v11.0.667 foi utilizado para a obtenção de um sumário mostrando a quantidade de leituras que passam em um filtros considerando limiares de estimativas de erros esperados (0,5, 1 e 2) em função de tamanhos de leituras (200pb a 300pb, com incremento de 10pb). Informações úteis para escolhas de parâmetros relacionados à poda de qualidade e fusão dos pares (R1 e R2) quando necessário. O módulo *fastq_eestats* do programa usearch v11.0.667 também foi utilizado para a obtenção de relatórios de qualidade com estimativas de erro esperado. A documentação das colunas pode ser consultada no [site do módulo](https://www.drive5.com/usearch/manual/cmd_fastq_eestats.html). ```{bash, label='usearch-eestats-reads', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/input" outdir="${BASE_DIR}/output/usearch/pre/trimsize" mkdir -p ${outdir} for fq in $(ls ${indir}/*.fastq) ; do bn=$(basename ${fq} .fastq) echo "Estimating errors according to trimming sizes in sample reads ${bn}" usearch11 \ -fastq_eestats2 ${fq} \ -output ${outdir}/${bn}_trimSizes.txt \ -length_cutoffs 200,300,10 \ -ee_cutoffs 0.5,1,2 > /dev/null 2> /dev/null done outdir="${BASE_DIR}/output/usearch/pre/eestats" mkdir -p ${outdir} for fq in $(ls ${indir}/*.fastq) ; do echo "Reporting expected error estimates in sample reads ${bn}" usearch11 \ -fastq_eestats ${fq} \ -output ${outdir}/${bn}_eestats.txt > /dev/null 2> /dev/null done ``` #### Resultados das estimativas de tamanhos e qualidades das leituras [Arquivos \*\_trimSizes.txt](./usearch/pre/trimsize/) [Arquivos \*\_eestats.txt](./usearch/pre/eestats/) ### Estimativas das médias de erros esperados (*reads*) Utilização do módulo *fastx_info* do programa usearch v11.0.667 foi utilizado para a obtenção da estimativa de erro médio dentro de cada biblioteca, considerando as suas *reads*. ```{bash, label='usearch-meanee-reads', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/input" outdir="${BASE_DIR}/output/usearch/pre/meanee" mkdir -p ${outdir} tablefile=${outdir}/meanEE.txt echo -e "ID\tEE" > ${tablefile} for fq in $(ls ${indir}/*.fastq) ; do bn=$(basename ${fq} .fastq) echo "Capturing 'Expected Error mean' of ${bn}" usearch11 \ -fastx_info ${fq} \ -output ${outdir}/${bn}_meanEE.txt > /dev/null 2> /dev/null ID=${bn} EE=$(cat ${outdir}/${bn}_meanEE.txt | grep -i "EE" | cut -d ' ' -f 3 | sed 's/;//g') echo -e "${ID}\t${EE}" >> ${tablefile} done ``` #### Resultado das estimativas das médias de erros esperados (*reads*) [Arquivo meanEE.txt](./usearch/pre/meanee/meanEE.txt) ### Podas de bases de baixa qualidade Em seguida, as leituras foram submetidas a podas de bases com baixa qualidade com o programa fastp v0.23.2, considerando a estratégia de janela deslizante, avaliando 5 bases a partir da extremidade direita, deslizando 1 base à esquerda e podando as bases da janela se a qualidade média dessas bases estiver abaixo de 10). O programa também foi configurado para descartar as leituras com a média de qualidade das bases abaixo de 10. As podas foram feitas independentemente para as leituras R1 e R2 de cada amostra. ```{bash, label='fastp-quality-trimming', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/input" outdir="${BASE_DIR}/output/fastp" mkdir -p ${outdir} for fq in $(ls ${indir}/*_R1.fastq) ; do bn=$(basename ${fq} _R1.fastq) echo "Quality evaluation of libraries ${bn}" fastp -i ${indir}/${bn}_R1.fastq \ -I ${indir}/${bn}_R2.fastq \ -o ${outdir}/${bn}_cleaned_R1.fastq \ -O ${outdir}/${bn}_cleaned_R2.fastq \ --disable_adapter_trimming \ --average_qual 10 \ --cut_right \ --cut_right_window_size 5 \ --cut_right_mean_quality 10 \ --html ${outdir}/${bn}_R1.html \ --json ${outdir}/${bn}_R1.json \ > ${outdir}/${bn}_cleaned.out.log \ 2>${outdir}/${bn}_cleaned.err.log done ``` ### Avaliação de qualidade das *reads* posterior à etapa de poda de qualidade A avaliação de qualidade das *reads* foi realizada da mesma forma como descrito anteriormente. ```{bash, label='fastqc-after', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/output/fastp" outdir="${BASE_DIR}/output/fastqc/fastp" mkdir -p ${outdir} rm -f ${outdir}/R1.txt rm -f ${outdir}/R2.txt for fq in $(ls ${indir}/*_R1*.fastq) ; do bn=$(basename ${fq} _R1.fastq) echo "Running FastQC for ${bn}_R1 and ${bn}_R2 ..." fastqc \ ${indir}/${bn}_R1.fastq \ ${indir}/${bn}_R2.fastq \ --threads 10 \ --outdir ${outdir} > /dev/null 2> /dev/null echo "${outdir}/${bn}_R1_fastqc.zip" >> ${outdir}/R1.txt echo "${outdir}/${bn}_R2_fastqc.zip" >> ${outdir}/R2.txt done multiqc --title "Post-processing reads (R1)" --file-list ${outdir}/R1.txt -o ${outdir} multiqc --title "Post-processing reads (R2)" --file-list ${outdir}/R2.txt -o ${outdir} ``` #### Resultados da avaliação de qualidade das *reads* após a poda de qualidade [MultiQC Report - R1](./fastqc/fastp/Post-processing-reads-R1_multiqc_report.html) [MultiQC Report - R2](./fastqc/fastp/Post-processing-reads-R2_multiqc_report.html) #### Fusão das *reads* R1 e R2 em sequências de *amplicons* As leituras R1 e R2 do sequenciamento foram fundidas para a representação dos *amplicons*. Para isso foi utilizado o programa flash v1.2.11, considerando mínimo de sobreposição de 15 bases (*--min-overlap 15*) e um máximo de 250 bp (*--max-overlap 250*), além de um máximo de falhas na correspondência de 20% (*--max-mismatch-density 0.2*). ```{bash, label='flash-merge', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/output/fastp" outdir="${BASE_DIR}/output/flash/pos" mkdir -p ${outdir} for fq in $(ls ${indir}/*_cleaned_R1.fastq) ; do bn=$(basename ${fq} _cleaned_R1.fastq) echo "Merge R1 and R2 reads of library ${bn}" flash \ ${indir}/${bn}_cleaned_R1.fastq \ ${indir}/${bn}_cleaned_R2.fastq \ --min-overlap 15 \ --max-overlap 250 \ --threads 10 \ --max-mismatch-density 0.2 \ --output-directory ${outdir} \ --output-prefix ${bn} > ${outdir}/${bn}.out.log 2> ${outdir}/${bn}.err.log done ``` ### Obtenção dos *amplicons* utilizando os *primers* especificados A obtenção dos *amplicons* via PCR *in-silico* foi realizada com o módulo *search_pcr* do programa usearch v11.0.667. Utilizando o *primer* *forward* e o *primer* *reverse*, que permitem a captura da região V3-V4 do gene 16S rRNA. Um máximo de 5 diferenças foi considerado como uma correspondência (*match*) aceita entre as sequências dos *primers* e as sequências dos *amplicons*. ```{bash, label='get-amplicons', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/output/flash/pos" outdir="${BASE_DIR}/output/usearch/pos/amplicons" mkdir -p ${outdir} minamp=300 maxamp=500 for fq in `ls ${indir}/*.extendedFrags.fastq 2>/dev/null`; do bn=`basename ${fq} .extendedFrags.fastq`; usearch11 -search_pcr ${fq} \ -db ${PRIMERS_FILE} \ -strand both \ -maxdiffs 5 \ -minamp ${minamp} \ -maxamp ${maxamp} \ -pcrout ${outdir}/${bn}.txt \ -ampoutq ${outdir}/${bn}.fastq done ``` Quando há mais do que um *amplicon* por leitura, então a maior é a selecionada com *in-house Perl script* [GetLargestAmp.pl](https://github.com/dgpinheiro/bioinfoutilities/blob/master/GetLargestAmp.pl). ```{bash, label='get-largest-amplicons', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}//output/usearch/pos/amplicons" outdir="${BASE_DIR}/output/amplicons" mkdir -p ${outdir} for fq in `ls ${indir}/*.fastq 2>/dev/null`; do bn=`basename ${fq} .fastq`; GetLargestAmp.pl -i ${fq} -o ${outdir}/${bn}.fastq done ``` ### Avaliação de qualidade dos *amplicons* obtidos após PCR *in-silico* A avaliação de qualidade dos *amplicons* foi realizada da mesma forma como descrito anteriormente para as *reads*. ```{bash, label='fastqc-amplicons', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/output/amplicons" outdir="${BASE_DIR}/output/fastqc/pos" mkdir -p ${outdir} for fq in $(ls ${indir}/*.fastq) ; do bn=$(basename ${fq} .fastq) echo "Running FastQC for ${bn} ..." fastqc \ ${indir}/${bn}.fastq \ --threads 15 \ --outdir ${outdir} > /dev/null 2> /dev/null done multiqc ${outdir} --title "Amplicons" -o ${outdir} ``` #### Resultados da avaliação de qualidade dos *amplicons* [MultiQC Report - Amplicons](./fastqc/pos/Amplicons_multiqc_report.html) ### Estimativas de tamanhos e qualidades das leituras (*amplicons*) Procedimento realizado exatamente como descrito anteriormente para as *reads*. ```{bash, label='usearch-eestats-amplicons', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/output/amplicons" outdir="${BASE_DIR}/output/usearch/pos/trimsize" mkdir -p ${outdir} for fq in $(ls ${indir}/*.fastq) ; do bn=$(basename ${fq} .fastq) echo "Estimating errors according to trimming sizes in sample amplicons ${bn}" usearch11 \ -fastq_eestats2 ${fq} \ -output ${outdir}/${bn}_trimSizes.txt \ -length_cutoffs 200,300,10 \ -ee_cutoffs 0.5,1,2 > /dev/null 2> /dev/null done outdir="${BASE_DIR}/output/usearch/pos/eestats" mkdir -p ${outdir} for fq in $(ls ${indir}/*.fastq) ; do echo "Reporting expected error estimates in sample amplicons ${bn}" usearch11 \ -fastq_eestats ${fq} \ -output ${outdir}/${bn}_eestats.txt > /dev/null 2> /dev/null done ``` #### Resultados das estimativas de tamanhos e qualidades das leituras [Arquivos \*\_trimSizes.txt](./usearch/pos/trimsize/) [Arquivos \*\_eestats.txt](./usearch/pos/eestats/) ### Estimativas das médias de erros esperados (*reads*) Utilização do módulo *fastx_info* do programa usearch v11.0.667 foi utilizado para a obtenção da estimativa de erro médio dentro de cada biblioteca, considerando as suas *reads*. ```{bash, label='usearch-meanee-amplicons', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/output/amplicons" outdir="${BASE_DIR}/output/usearch/pos/meanee" mkdir -p ${outdir} tablefile=${outdir}/meanEE.txt echo -e "ID\tEE" > ${tablefile} for fq in $(ls ${indir}/*.fastq) ; do bn=$(basename ${fq} .fastq) echo "Capturing 'Expected Error mean' of ${bn}" usearch11 \ -fastx_info ${fq} \ -output ${outdir}/${bn}_meanEE.txt > /dev/null 2> /dev/null ID=${bn} EE=$(cat ${outdir}/${bn}_meanEE.txt | grep -i "EE" | cut -d ' ' -f 3 | sed 's/;//g') echo -e "${ID}\t${EE}" >> ${tablefile} done ``` #### Resultado das estimativas das médias de erros esperados (*amplicons*) [Arquivo meanEE.txt](./usearch/pos/meanee/meanEE.txt) ### Contagens preliminares ```{bash, label='fastq-counts', engine.opts='-l'} #| eval: false #| echo: false indir="${BASE_DIR}/input" outdir="${BASE_DIR}/output" mkdir -p "${outdir}/info" table=${outdir}/info/stats.txt rm -f ${table} echo -e "ID\tRaw\tCleaned\tMerged\tValid" > ${table} for fastq in $(ls ${indir}/*_R1.fastq) ; do bn=$(basename ${fastq} _R1.fastq) id=$(echo -e ${bn}) raw=$(cat ${fastq} | wc -l | sed 's/ /\t/g' | awk '{print $1/4}') cln=$(cat ${outdir}/fastp/${bn}_cleaned_R1.fastq | wc -l | sed 's/ /\t/g' | awk '{print $1/4}') mgd=$(cat ${outdir}/flash/pos/${bn}.extendedFrags.fastq | wc -l | sed 's/ /\t/g' | awk '{print $1/4}') vld=$(cat ${outdir}/amplicons/${bn}.fastq | wc -l | sed 's/ /\t/g' | awk '{print $1/4}') echo -e "${id}\t${raw}\t${cln}\t${mgd}\t${vld}" >> ${table} done ``` ### Processamento dos Amplicons e obtenção das *Amplicon Sequence Variants* (ASVs) As etapas de processamento a seguir foram realizadas no ambiente R v4.1.2. ```{r, label="setup-dada2-main-path"} #| echo: false #| eval: true caminho <- paste0(base.dir, "/output/dada2/") ``` ```{r, label="setup-dada2-other-paths"} #| echo: false #| eval: false # Criando um atalho para o caminho que iremos salvar as saídas: if (! dir.exists(caminho) ) { dir.create(caminho) } input <- data.frame(reads=list.files(paste0(caminho,"../amplicons"), pattern = ".fastq", full.names=TRUE)) input$nomes <- sapply(strsplit(basename(input$reads), "\\."), `[`, 1) rownames(input) <- input$nomes # Separando nome das amostras: nomes <- input$nomes # Carregando os caminhos dos arquivos reads <- input$reads names(reads) <- nomes # Criando caminhos de saída para as bibliotecas processadas com dada2: reads_filt <- file.path(caminho, "filterAndTrim", paste0(nomes, ".fastq")) names(reads_filt) <- nomes ``` ```{r, label="load-R-environment"} #| echo: false #| include: false #| warning: false #| eval: true options(getClass.msg=FALSE) rand.samp <- NULL load(file = paste0(caminho,'/Analysis_Amplicon_RDP.RData')) if (is.null(rand.samp)) { rand.samp <- sample(1:length(reads),1) } ``` #### Avaliação do perfil de qualidade das sequências dos amplicons As leituras R1 foram carregadas no pacote DADA2 para o reconhecimento dos amplicons. Na @fig-plotQPr a seguir temos uma análise da qualidade das leituras em função do tamanho. ```{r, label="plotQualityProfile-reads"} #| echo: false #| warning: false #| eval: false # QC plotQPr <- suppressWarnings( plotQualityProfile(reads) ) # Visualização pré-QC svg(filename = paste0(caminho,'Amplicons_Quality_Score_pre_DADA2','.svg'), width = 20, height = 20, pointsize = 12 ) plotQPr graphics.off() plotQPr_sample <- suppressWarnings( plotQualityProfile(reads[rand.samp]) ) # Visualização pré-QC svg(filename = paste0(caminho,'Amplicons_Quality_Score_pre_DADA2_sample','.svg'), width = 20, height = 20, pointsize = 12 ) plotQPr_sample graphics.off() ``` ```{r, label="plotQualityProfile-reads-shown"} #| echo: false #| warning: false #| lightbox: true #| fig.align: center #| out.width: 50% #| label: fig-plotQPr #| fig-cap: "Quality Profile previous DADA2 processing (only one SAMPLE)" #| eval: true knitr::include_graphics(paste0(caminho,'Amplicons_Quality_Score_pre_DADA2_sample','.svg')) ``` [Quality Profile previous DADA2 processing (full image)](./dada2/Amplicons_Quality_Score_pre_DADA2.svg) #### Filtragem e poda de sequências de amplicons A seguir foi realizada uma filtragem e poda nas sequências de acordo com o perfil de qualidade observado, dessa forma, consideramos não truncar as sequências até um determinado tamanho (*truncLen = 0*), com um valor máximo de expectativa de erros (*maxEE*) de 0,5, removendo amplicons provenientes do controle PhiX e as de baixa complexidade (número efetivo de *kmers* menor que 5). A intenção também é a de manter o máximo de leituras possível, sem a perda de qualidade. ```{r, label="filterAndTrim"} #| echo: false #| eval: false #TO_UNCOMMENT qc <- filterAndTrim(fwd = reads, filt = reads_filt, truncLen = 180, trimLeft = 0, maxEE = 0.5, maxN = 0, rm.phix = T, compress = F, multithread=T, rm.lowcomplex=5, verbose=T ) rownames(qc) <- gsub('.fastq','',rownames(qc)) ``` O resultado (@tbl-filttrim) pode também ser visualizado na @fig-plotQPrf a seguir para apenas uma amostra `r names(reads[rand.samp])`, como representante das demais. O download da imagem completa com todas as amostras pode ser feito pelos links seguintes às imagens. ```{r, label="filterAndTrim-shown"} #| echo: false #| label: tbl-filttrim #| tbl-cap: Read after filtering and trimming process #| eval: true knitr::kable(qc) ``` ```{r, label="plotQualityProfile-filtered-reads"} #| echo: false #| warning: false #| eval: false plotQPrf <- suppressWarnings( plotQualityProfile(reads_filt) ) # Visualização pós-QC svg(filename = paste0(caminho,'Amplicons_Quality_Score_DADA2_filterAndTrim','.svg'), width = 20, height = 20, pointsize = 12 ) plotQPrf graphics.off() plotQPrf_sample <- suppressWarnings( plotQualityProfile(reads_filt[rand.samp]) ) # Visualização pós-QC svg(filename = paste0(caminho,'Amplicons_Quality_Score_DADA2_filterAndTrim_sample','.svg'), width = 20, height = 20, pointsize = 12 ) plotQPrf_sample graphics.off() ``` ```{r, label="plotQualityProfile-filtered-reads-shown"} #| echo: false #| warning: false #| lightbox: true #| fig.align: center #| out.width: 50% #| label: fig-plotQPrf #| fig-cap: "Quality Profile post DADA2 filtering and trimming (only one sample)" #| eval: true knitr::include_graphics(paste0(caminho,'Amplicons_Quality_Score_DADA2_filterAndTrim_sample','.svg')) ``` [Quality Profile post DADA2 filtering and trimming (full image)](Amplicons_Quality_Score_DADA2_filterAndTrim.svg) #### Identificação das taxas de erros Os amplicons passaram então por um processo de identificação das taxas de erros a partir de um processo de aprendizado (*learnErrors*) que alterna entre estimativa de erros e inferência de um modelo para os erros apartir dos dados até atingir convergência. A representação gráfica do resultado deste processo pode ser observado na @fig-plotErrors. A sugestão proposta pelo autor do dada2 de forçar a monotonicidade do modelo de erros ajustado foi seguida. O número máximo de iterações na etapa de autoconsistência para estimar a taxa de erros foi incrementada para 25 (MAX_CONSIST = 25). A convergência foi atingida após 15 iterações. ```{r, label="learnErrors"} #| echo: false #| warning: false #| eval: false erros <- learnErrors(reads_filt, multithread = T, randomize=T, verbose=T, MAX_CONSIST = 25, errorEstimationFunction=loessErrfun) ``` ```{r, label="plotErrors"} #| echo: false #| include: false #| eval: false plotErr1 <- plotErrors(erros, nominalQ = TRUE) ggsave(paste0(caminho,'Amplicons_Quality_Score_DADA2_learnErrors_1','.svg'), plotErr1, width=12, unit='cm' ) erros.md.full <- erros ## sugerido pelo autor do dada2 ## https://github.com/benjjneb/dada2/issues/791 make.monotone.decreasing <- function(v) sapply(seq_along(v), function(i) max(v[i:length(v)])) erros.md <- t(apply(getErrors(erros), 1, make.monotone.decreasing)) erros.md.full <- erros colnames(erros.md) <-colnames(erros.md.full$err_out) erros.md.full$err_out <- erros.md plotErr2 <- plotErrors(erros.md.full, nominalQ = TRUE) ggsave(paste0(caminho,'Amplicons_Quality_Score_DADA2_learnErrors_2','.svg'), plotErr2, width=12, unit='cm' ) ``` ```{r, label="plotErrors-shown"} #| echo: false #| warning: false #| lightbox: true #| fig.align: center #| out.width: 50% #| label: fig-plotErrors #| fig-cap: "The error rates for each possible transition (A→C, A→G, …) are shown. Points are the observed error rates for each consensus quality score. The black line shows the estimated error rates after convergence of the machine-learning algorithm. The red line shows the error rates expected under the nominal definition of the Q-score. " #| eval: true knitr::include_graphics(paste0(caminho,'Amplicons_Quality_Score_DADA2_learnErrors_2','.svg')) ``` #### Inferência das *Amplicon Sequence Variants* (ASVs) As sequências de amplicons devidamente podadas e filtradas foram submetidas à deduplicação (*dereplication*) para em seguida o algoritmo do DADA2 inferir as *Amplicon Sequence Variants* (ASVs) seguindo um processo de remoção de ruídos (*denoising*) que leva em conta as taxas de erros identificadas na etapa anterior. ```{r, label="dada2"} #| echo: false #| eval: false reads_derep <- derepFastq(reads_filt) # Deduplicação reads_dada <- dada(reads_derep, err = erros.md.full, multithread = T, pool=FALSE) # Denoising rm(reads_derep) asvs <- makeSequenceTable(reads_dada) ``` #### Identificação e filtragem de ASVs quiméricas O resultado foi de `r dim(asvs)[2]` ASVs. As quais foram submetidas à identificação de quimeras (*removeBimeraDenovo*) pelo método *consensus* (as sequências das amostras são analisadas independentemente e uma decisão é tomada por consenso para cada ASV), considerando um mínimo de 100 de abundância para as ASVs que podem ser consideradas candidatas a aparentadas (origem das quimeras) além de uma razão de abundância de no mínimo 10 em relação às potenciais quimeras. ```{r, label="removeBimeraDenovo"} #| echo: false #| include: false #| eval: false asvs_nochim <- removeBimeraDenovo(asvs, method = "consensus", multithread = T, minParentAbundance=100, allowOneOff=T, minFoldParentOverAbundance=10, verbose=T) ``` Foram identificadas `r dim(asvs)[2]-dim(asvs_nochim)[2]` ASVs quiméricas. #### Números de leituras iniciais e abundâncias de ASVs Abaixo (@tbl-readsasvs) os números de leituras iniciais e de abundâncias de ASVs após os processamentos até o momento. ```{r, label="preliminary-stats"} #| echo: false #| include: false #| eval: false # Resumo dos processamentos no DADA cont_dada <- data.frame(cbind(nomes, qc[,2], sapply(reads_dada, getN), rowSums(asvs_nochim))) colnames(cont_dada) <- c("ID", "filterAndTrim", "Denoised", "ASVs") rownames(cont_dada) <- NULL # Carregando resumo do pre-processamento: cont_preproc <- read.table(paste0(caminho,"../info/stats.txt"), header = T, sep = "\t") # Fundindo os resumos preproc + DADA: cont_proc <- merge(cont_preproc, cont_dada, by = "ID") %>% mutate_at(-1, as.numeric) cont_proc.plot <- select(cont_proc, -c(Valid)) colnames(cont_proc.plot) <- c('ID','Raw.reads','Cleaned.reads','FilteredAndTrimmed.reads','Denoised.asvs.abundance','ChimeraFree.asvs.abundance') # Salvando tabela write.table(cont_proc.plot, file = paste0(caminho,"../info/summary.txt"), sep = "\t", col.names = T, row.names = F) ``` ```{r, label="preliminary-stats-shown"} #| echo: false #| label: tbl-readsasvs #| tbl-cap: Reads and ASV abundances #| eval: true knitr::kable(cont_proc.plot) ``` #### Tabela de contagem ASVs ```{r, label="ASVs-fasta"} #| echo: false #| eval: false # Obtendo tabela de contagem: contagens_df <- data.frame(t(asvs_nochim)) %>% mutate(seq = rownames(.)) %>% select(seq, everything()) writeFasta(contagens_df[,'seq'], file=paste0(caminho,"ASVs.fa")) ``` O número total de ASVs válidas é de `r dim(contagens_df)[1]`. As sequências fasta das ASVs podem ser acessadas no link abaixo. ```{r, label="ASVs-fasta-shown"} #| echo: false #| eval: true xfun::embed_file(paste0(caminho,'ASVs','.fa')) ``` #### Atribuição taxonômica às ASVs Para a atribuição taxonômica, as funções *assignTaxonomy()* e *addSpecies()*, ambas do pacote DADA2 foram utilizadas com base nos bancos de dados GTDB (*GTDB_bac120_arc122_ssu_r202*), RDP (*rdp_train_set_18*), RefRDP (*RefSeq_16S_6-11-20_RDPv16*), e SILVA (*silva_nr99_v138.1_train_set*). As possibilidades de correspondências foram consideradas em ambas as orientações e com um mínimo de valor de bootstrap de 60%. ```{r, label="taxonomy-assignment"} #| echo: false #| eval: false ## GTDB (Genome Taxonomy Database): taxGTDB <- assignTaxonomy(asvs_nochim, "/data/DB/taxonomy/GTDB_bac120_arc122_ssu_r202_fullTaxo.fa.gz", multithread = T, tryRC = T, minBoot=80) taxGTDB <- addSpecies(taxGTDB, "/data/DB/taxonomy/GTDB_bac120_arc122_ssu_r202_Species.fa.gz", tryRC = T) ## RDP (Ribosomal Database Project): taxRDP <- assignTaxonomy(asvs_nochim, "/data/DB/taxonomy/rdp_train_set_18.fa.gz", multithread = T, tryRC = T, minBoot=80) taxRDP <- addSpecies(taxRDP, "/data/DB/taxonomy/rdp_species_assignment_18.fa.gz", tryRC = T) ## RefRDP (NCBI RefSeq 16S rRNA database supplemented by RDP): taxRefRDP <- assignTaxonomy(asvs_nochim, "/data/DB/taxonomy/RefSeq_16S_6-11-20_RDPv16_fullTaxo.fa.gz", multithread = T, tryRC = T, minBoot=80) taxRefRDP <- addSpecies(taxRefRDP, "/data/DB/taxonomy/RefSeq_16S_6-11-20_RDPv16_Species.fa.gz", tryRC = T) ## SILVA: taxSilva <- assignTaxonomy(asvs_nochim, "/data/DB/taxonomy/silva_nr99_v138.1_train_set.fa.gz", multithread = T, tryRC = T, minBoot=80) taxSilva <- addSpecies(taxSilva, "/data/DB/taxonomy/silva_species_assignment_v138.1.fa.gz", tryRC = T) ## Transformações das tabelas de taxonomia: taxGTDB_df <- taxGTDB %>% data.frame() %>% # Transformando de matriz para data frame mutate(seq = rownames(.)) %>% # Criando uma coluna para a ASV select(seq, Kingdom:Species) %>% # Alterando a ordem das colunas mutate(Species = str_remove_all(Species, "\\(.*$")) %>% # Removendo alguns apêndices das espécies. Ids, cepas... mutate(Species = str_replace_all(Species, " ", "_")) %>% # Trocando espaços por "_" mutate(across(Kingdom:Species, ~ gsub("^$", NA, .))) # Trocando células vazias por "NA" taxRDP_df <- taxRDP %>% data.frame() %>% mutate(seq = rownames(.)) %>% mutate(Species2 = paste0(Genus, "_", Species)) %>% mutate(Species2 = str_remove_all(Species2, ".*NA$")) %>% mutate(Species = Species2) %>% select(seq, Kingdom:Species) %>% mutate(across(Kingdom:Species, ~ gsub("^$", NA, .))) taxRefRDP_df <- taxRefRDP %>% data.frame() %>% mutate(seq = rownames(.)) %>% select(seq, Kingdom:Species) %>% mutate(across(Kingdom:Species, ~ gsub("^$", NA, .))) %>% mutate(Species = str_remove_all(Species, "\\(.*$")) %>% mutate(Species = str_replace_all(Species, " ", "_")) taxSilva_df <- taxSilva %>% data.frame() %>% mutate(seq = rownames(.)) %>% mutate(Species2 = paste0(Genus, "_", Species)) %>% mutate(Species2 = str_remove_all(Species2, ".*NA$")) %>% mutate(Species = Species2) %>% select(seq, Kingdom:Species) %>% mutate(across(Kingdom:Species, ~ gsub("^$", NA, .))) ##Tabela de contagens completa taxAll_cont_df <- merge(taxGTDB_df %>% mutate(GTDB = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, GTDB), merge(taxRDP_df %>% mutate(RDP = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, RDP), merge(taxRefRDP_df %>% mutate(RefRDP = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, RefRDP), merge(taxSilva_df %>% mutate(Silva = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, Silva), contagens_df, by = "seq"), by = "seq"), by = "seq"), by = "seq") ## Estabelecendo contaminantes: contaminantes <- c("Chloroplast;", "Mitochondria;", "^Unclassified;", "^Eukaryota;") ## Extraindo seqs. não contaminantes: seqs_boas <- taxAll_cont_df %>% filter_at(vars(-seq), all_vars(str_detect(., paste(contaminantes, collapse = "|"), negate = T))) %>% pull(seq) ## Filtrando seqs. contaminantes: taxGTDB_filt_df <- taxGTDB_df %>% filter(seq %in% seqs_boas) taxRDP_filt_df <- taxRDP_df %>% filter(seq %in% seqs_boas) taxRefRDP_filt_df <- taxRefRDP_df %>% filter(seq %in% seqs_boas) taxSilva_filt_df <- taxSilva_df %>% filter(seq %in% seqs_boas) contagens_filt_df <- contagens_df %>% filter(seq %in% seqs_boas) ## Comparação das classificações - Função customizada: ClassTax <- summary.tax(ordered_names = c("GTDB", "RDP", "RefRDP", "Silva"), taxGTDB_filt_df, taxRDP_filt_df, taxRefRDP_filt_df, taxSilva_filt_df) sample_cols <- setdiff(colnames(contagens_filt_df), 'seq') ## Seqs. classificadas até determinado nível: ClassTax_cont <- ClassTax[[1]] ## Táxons únicos até determinado nível: ClassTax_unic <- ClassTax[[2]] ``` ##### Identificação e filtragem de ASVs contaminantes Além disso foram feitas filtragens de ASVs contaminantes ou não classificadas (Padrões buscados: `r str_escape( paste(gsub('(.*)','"\\1"',as.character(contaminantes)),collapse=","))`) resultando em `r length(seqs_boas)` ASVs válidas livres de contaminantes. ```{r, label="ContaminantFree_ASVs-fasta"} #| echo: false #| eval: false writeFasta(seqs_boas, file=paste0(caminho,'ContaminantFree_ASVs','.fa')) ``` As sequências fasta livre de contaminantes podem ser acessadas no link abaixo. ```{r, label="ContaminantFree_ASVs-fasta-shown"} #| echo: false #| eval: true xfun::embed_file(paste0(caminho,'ContaminantFree_ASVs','.fa')) ``` ##### Resultados das atribuições taxonômicas As @tbl-classtax_cont e @tbl-classtax_unic sumarizam o resultado das atribuições taxonômicas correspondentes aos diferentes bancos de dados considerados no processo. ```{r, label="ClassTax_cont-shown"} #| echo: false #| label: tbl-classtax_cont #| tbl-cap: Somatório da abundância de ASVs com atribuição até determinado nível taxonômico #| eval: true knitr::kable(ClassTax_cont) ``` ```{r, label="ClassTax_uniq-shown"} #| echo: false #| label: tbl-classtax_unic #| tbl-cap: ASVs com atribuição até determinado nível taxonômico #| eval: true knitr::kable(ClassTax_unic) ``` ##### Planilha das ASVs com atribuição taxonômica {#sec-planilha} Por fim, selecionamos as atribuições taxonômicas realizadas com o banco de dados RDP, afinal é um banco de dados que possui curadoria, e portanto, confiabilidade. Além disson, possui uma quantidade satisfatória de ASVs reconhecidas até o nível de gênero. A planilha contendo os resultados da atribuição taxonômica com o banco de dados RDP pode ser acessada no link abaixo. Além das atribuições taxonômicas ao longo dos níveis, a planilha contém a soma da abundância e o número que representa a prevalência das ASVs entre as amostras. ```{r, label="master-table-write"} #| echo: false #| eval: false ## Construção da "master table": master <- merge(taxRDP_filt_df, contagens_filt_df) %>% # Fusão das taxonomias e contagens mutate("CountSum" = rowSums(.[9:ncol(.)])) %>% # Somatória de abundancia de cada ASV mutate("CountPrev" = rowSums(.[9:ncol(.)] != 0) - 1) %>% # Somatória de prev. nas amostras (descontando a somatória) arrange(-CountSum) %>% # Ordenando tabela pela abundancia mutate("ID" = paste0("ASV_", seq(1, nrow(.)))) %>% # Adicionando identificador único select(ID, everything()) # Reordenando... file.remove(paste0(caminho,'master.xlsx')) XLConnect::xlcFreeMemory() options(java.parameters = "-Xmx8192m") openxlsx::write.xlsx(x = master, file=paste0(caminho,'master.xlsx'), sheetName = "Master") ``` ```{r, label="master-table-shown"} #| echo: false #| eval: true xfun::embed_file(paste0(caminho,'master.xlsx')) ``` ```{r, label="taxonomy-assignments-percents"} #| echo: false #| eval: false # Porcentagens de classificações: fam_perc <- sum(master %>% filter(!is.na(Family)) %>% pull(CountSum)) / sum(master$CountSum)*100 # % fam gen_perc <- sum(master %>% filter(!is.na(Genus)) %>% pull(CountSum)) / sum(master$CountSum)*100 # % gen spc_perc <- sum(master %>% filter(!is.na(Species)) %>% filter(Species!="") %>% pull(CountSum)) / sum(master$CountSum)*100 # % sp ``` No geral, houve `r round(fam_perc,2)`% de dados com identificação da família, `r round(gen_perc,2)`% de dados com informação de gênero, e `r round(spc_perc,2)`% de dados com informação de espécie. A tabela anterior (@tbl-readsasvs) com os números de leituras portanto foi atualizada (@tbl-readsasvs-up) com novas colunas. ```{r, label="summary-write"} #| echo: false #| include: false #| eval: false ## Contagens finais: cont_proc$ASVs_clean <- master %>% select(all_of(sample_cols)) %>% colSums() cont_proc$Usable_perc <- round((cont_proc$ASVs_clean / cont_proc$Raw) * 100, 2) cont_proc.plot$ContaminantFree.asvs.abundance <- cont_proc$ASVs_clean cont_proc.plot$PercentUsable.reads <- cont_proc$Usable_perc if (!dir.exists(paste0(caminho,'info'))) { dir.create(paste0(caminho,'info')) } write.table(cont_proc.plot, file = paste0(caminho,'info/summary.txt'), sep = "\t", col.names = T, row.names = F) ``` ```{r, label="summary-shown"} #| echo: false #| label: tbl-readsasvs-up #| tbl-cap: Reads and ASVs abundances (updated) #| eval: true knitr::kable(cont_proc.plot) ``` ##### Análise Filogenética As `r length(seqs_boas)` livres de contaminantes foram alinhadas com pacote *DECIPHER* e usando pacote *phangorn* foi estimada uma árvore filogenética de Máxima Verossimilhança com a função *pml()* , a partir de uma árvore Neighbor-Joining utilizando parâmetros pré-definidos (*default*). Essa árvore inicial foi otimizada com a função *optimim.pml()* considerando rearranjos nos ramos por NNI (*Nearest Neighbor Interchange*), o modelo evolutivo General Time Reversible (GTR), e a otimização de parâmetros do modelo, tais como a proporção de sítios invariáveis (a partir do valor 0,2) e o parâmetro gamma. Em resumo, a análise filogenética segue a sugestão de [@callahan2016] exceto pelo método de rearranjo dos ramos. ```{r, label="phylogenetic-analysis"} #| echo: false #| warnings: false #| include: false #| eval: false ## https://compbiocore.github.io/metagenomics-workshop/assets/DADA2_tutorial.html ## https://www.mimuw.edu.pl/~lukaskoz/teaching/sad2/books/Analysis_of_Phylogenetics_and_Evolution_with_R.pdf ## https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4955027/ #Run Sequence Alignment (MSA) using DECIPHER seqs_boas <- master$seq names(seqs_boas) <- master$ID alignment <- AlignSeqs(DNAStringSet(seqs_boas), anchor=NA) #Change sequence alignment output into a phyDat structure phang.align <- phyDat(as(alignment, "matrix"), type="DNA") #Create distance matrix dm <- dist.ml(phang.align) #Perform Neighbor joining treeNJ <- NJ(dm) # Note, tip order != sequence order #Internal maximum likelihood fit = pml(treeNJ, data=phang.align) #negative edges length changed to 0! fitGTR <- update(fit, k=4, inv=0.2) fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE, rearrangement = "NNI", control = pml.control(trace = 0)) #rm(alignment) #rm(phang.align) #rm(dm) #rm(fit) #rm(treeNJ) ``` #### Arquivos para exploração com MicrobiomeAnalyst Os arquivos a seguir servem para análises em plataformas externas como o [MicrobiomeAnalyst](https://www.microbiomeanalyst.ca/MicrobiomeAnalyst/upload/OtuUploadView.xhtml "MicrobiomeAnalyst - Marker Data Profiling"). Essa plataforma exige que os dados estejam em replicatas (NO MÍNIMO 3 RÉPLICAS), dessa forma, cada amostra será replicada 3 vezes por um procedimento de subamostragem com repetição. O procedimento de subamostragem foi realizado com a função *rarefy_even_depth()* do pacote phyloseq v1.38.0 do R. ```{r, label="MicrobiomeAnalyst-files"} #| echo: false #| warnings: false #| include: false #| message: false #| eval: false # # # Salvando tabelas para o MicrobiomeAnalyst: tb_meta <- read.delim(file = paste0(base.dir,"/METADADOS_MicrobiomaFrangos.csv"), header = T, sep = "\t", check.names=FALSE) %>% dplyr::rename("#NAME"="ID") tb_cont <- master %>% select(all_of(c("ID", sample_cols))) %>% dplyr::rename("#NAME"="ID") tb_taxa <- master %>% select(ID, Kingdom:Species) %>% fix_tax() %>% # Função customizada dplyr::rename("#TAXONOMY"="ID") #View(tb_meta) #View(tb_cont) #View(tb_taxa) tb_cont.ma <- tb_cont %>% rename_at( vars(nomes), ~ gsub("(.*)","\\1.1",nomes) ) %>% bind_cols(tb_cont[,-1]) %>% rename_at( vars(nomes), ~ gsub("(.*)","\\1.2",nomes) ) %>% bind_cols(tb_cont[,-1]) %>% rename_at( vars(nomes), ~ gsub("(.*)","\\1.3",nomes) ) %>% remove_rownames() %>% tibble::column_to_rownames(var = "#NAME") tb_meta %>% mutate(`#NAME`= gsub("(.*)","\\1.1",`#NAME`) ) %>% bind_rows(tb_meta %>% mutate(`#NAME`= gsub("(.*)","\\1.2",`#NAME`) ) ) %>% bind_rows(tb_meta %>% mutate(`#NAME`= gsub("(.*)","\\1.3",`#NAME`) ) ) %>% tibble::column_to_rownames("#NAME") -> tb_meta.ma taxonomy <- tb_taxa %>% tibble::column_to_rownames("#TAXONOMY") otu <- tb_cont %>% remove_rownames() %>% tibble::column_to_rownames(var = "#NAME") metadata <- tb_meta %>% tibble::column_to_rownames("#NAME") options(getClass.msg=FALSE) for (n in nomes) { project.ma <- phyloseq::phyloseq(phyloseq::otu_table(as.matrix(tb_cont.ma[,c(paste0(n,'.1'),paste0(n,'.2'),paste0(n,'.3'))]), taxa_are_rows=TRUE), phyloseq::tax_table(as.matrix(taxonomy)), phyloseq::sample_data(tb_meta.ma[c(paste0(n,'.1'),paste0(n,'.2'),paste0(n,'.3')),]), phyloseq::phy_tree(fitGTR$tree) ) set.seed(1) phy_tree(project.ma) <- root(phy_tree(project.ma), sample(taxa_names(project.ma), 1), resolve.root = TRUE) # rarefy with replacement project.rarefied.ma = suppressMessages( rarefy_even_depth( project.ma, rngseed=1, replace=T) ) for (asv in rownames( as.data.frame(otu_table(project.rarefied.ma)) ) ) { tb_cont.ma[asv, paste0(n,'.1')] <- as.data.frame(otu_table(project.rarefied.ma))[asv, paste0(n,'.1')] tb_cont.ma[asv, paste0(n,'.2')] <- as.data.frame(otu_table(project.rarefied.ma))[asv, paste0(n,'.2')] tb_cont.ma[asv, paste0(n,'.3')] <- as.data.frame(otu_table(project.rarefied.ma))[asv, paste0(n,'.3')] } rm(project.ma) rm(project.rarefied.ma) } tb_cont.ma <- tb_cont.ma %>% rownames_to_column(var = "#NAME") tb_meta.ma <- tb_meta.ma %>% rownames_to_column(var = "#NAME") write.table(tb_cont.ma, file = paste0(caminho, "MA_counts.txt"), quote = F, sep = "\t", col.names = T, row.names = F) write.table(tb_meta.ma, file = paste0(caminho, "MA_metadata.txt"), quote = F, sep = "\t", col.names = T, row.names = F) write.table(tb_taxa, file = paste0(caminho, "MA_taxonomy.txt"), quote = F, sep = "\t", col.names = T, row.names = F) write.tree(fitGTR$tree, file = paste0(caminho, "MA_tree.nwk")) ``` ```{r, label="MA_counts.txt"} #| echo: false #| warnings: false #| eval: true xfun::embed_file(paste0(caminho, "MA_counts.txt")) ``` ```{r, label="MA_metadata.txt"} #| echo: false #| warnings: false #| eval: true xfun::embed_file(paste0(caminho, "MA_metadata.txt")) ``` ```{r, label="MA_taxonomy.txt"} #| echo: false #| warnings: false #| eval: true xfun::embed_file(paste0(caminho, "MA_taxonomy.txt")) ``` ```{r, label="MA_tree.nwk"} #| echo: false #| warnings: false #| eval: true xfun::embed_file(paste0(caminho, "MA_tree.nwk")) ``` ### Análises ecológicas ```{r, label="creating-phyloseq-object"} #| echo: false #| warnings: false #| include: false #| message: false #| eval: false project.ps <- phyloseq::phyloseq(phyloseq::otu_table(as.matrix(otu),taxa_are_rows=TRUE), phyloseq::tax_table(as.matrix(taxonomy)), phyloseq::sample_data(subset(metadata)), phyloseq::phy_tree(fitGTR$tree) ) set.seed(1) phy_tree(project.ps) <- root(phy_tree(project.ps), sample(taxa_names(project.ps), 1), resolve.root = TRUE) #is.rooted(phy_tree(project.ps)) rarefy.threshold <- 1*as.integer(quantile(sample_sums(project.ps),0.25)) project.small.ps <- phyloseq::phyloseq(phyloseq::otu_table(as.matrix(otu),taxa_are_rows=TRUE), phyloseq::tax_table(as.matrix(taxonomy)), phyloseq::sample_data( metadata[names( sample_sums( project.ps )[sample_sums(project.ps) < rarefy.threshold] ), ] ), phy_tree(project.ps) ) # rarefy without replacement project.final.ps <- rarefy_even_depth( project.ps, sample.size=rarefy.threshold, replace=F, trimOTUs=F, rngseed=1, verbose=TRUE) project.final.ps <- merge_phyloseq(project.final.ps, project.small.ps ) rmtaxa <- phyloseq::taxa_names(project.final.ps)[phyloseq::taxa_sums( project.final.ps ) <= 0] if (length(rmtaxa) > 0) { project.final.ps <- prune_taxa(setdiff(phyloseq::taxa_names(project.final.ps), rmtaxa), project.final.ps) } phy_tree(project.final.ps) <- root(phy_tree(project.final.ps), sample(taxa_names(project.final.ps), 1), resolve.root = TRUE) summ.ps<-summarize_phyloseq(project.ps) summ.rarefied.ps<-summarize_phyloseq(project.final.ps) summ.final.ps<-summarize_phyloseq(project.final.ps) #sparsity #(sum(otu_table(project.final.ps)==0)/(dim(otu_table(project.final.ps))[1]*dim(otu_table(project.final.ps))[2]))*100 ``` ```{r, label="sparsity"} #| echo: false #| warnings: false #| include: false #| eval: true spa<- suppressMessages( gsub('^7] ','',summarize_phyloseq(project.ps)[[6]]) ) spa.rarefied <- suppressMessages( gsub('^7] ','',summarize_phyloseq(project.final.ps)[[6]]) ) ``` Utilizamos os pacotes phyloseq v1.38.0 e vegan v2.6-4 no R para as análises subsequentes. A seguir temos um resultado do quanto a matriz de ASVs é esparsa [^1] usando a função *summarize_phyloseq()* que sumariza os resultados do objeto phyloseq sem a rarefação (`r spa`). [^1]: esparsa: valores espalhados, ou seja, com muitos zeros na matriz de abundâncias. #### Rarefação dos dados Nessas análises consideramos os dados após a rarefação em relação ao 25º percentil da distribuição de tamanho das amostras (`r options(scipen = -0, digits = 4); as.character(as.integer(quantile(sample_sums(project.ps),0.25)))`) (@fig-rarecurve). A rarefação foi feita para as amostras com um número de contagens maiores que o valor definido para a rarefação. As amostras menores que o valor definido foram apenas acrescentadas. Por fim, as ASVs que após a rarefação ficaram zeradas (`r options(scipen = -0, digits = 4); as.character(as.integer( length(rmtaxa) ))`) foram removidas. Também podem ser consultados mais detalhes na análise de rarefação por grupos de amostras (@fig-rarecurve-groups). As figuras foram geradas a partir da suavização dos valores de contagem de *reads* e respectivos intervalos de confiança (95%) quando apresentados. ```{r, label="rarefaction-curves"} #| echo: false #| warnings: false #| include: false #| eval: false sex.item <- c("FEMALE","MALE") sex.item.color <- c("deeppink","steelblue") sex.item.linetype <- c("dashed","solid") treatment.item <- c("CONTROL","VEHICLE","PROBIOTIC","ANTIBIOTIC") treatment.item.color <- c("#BDC3C7","#85929E","#424949","#232429") treatment.item.linetype <- c('solid','solid','solid','solid') segment.item <- c("DUODENUM","JEJUNUM","ILEUM","CECUM") segment.item.color <- c("#1B4F72","#4A235A","#186A3B","#6E2C00") segment.item.linetype <- c('solid','solid','solid','solid') sample.item <- c( 'DUODENUM_CONTROL_FEMALE', 'DUODENUM_CONTROL_MALE', 'DUODENUM_VEHICLE_FEMALE', 'DUODENUM_VEHICLE_MALE', 'DUODENUM_PROBIOTIC_FEMALE', 'DUODENUM_PROBIOTIC_MALE', 'DUODENUM_ANTIBIOTIC_FEMALE', 'DUODENUM_ANTIBIOTIC_MALE', 'JEJUNUM_CONTROL_FEMALE', 'JEJUNUM_CONTROL_MALE', 'JEJUNUM_VEHICLE_FEMALE', 'JEJUNUM_VEHICLE_MALE', 'JEJUNUM_PROBIOTIC_FEMALE', 'JEJUNUM_PROBIOTIC_MALE', 'JEJUNUM_ANTIBIOTIC_FEMALE', 'JEJUNUM_ANTIBIOTIC_MALE', 'ILEUM_CONTROL_FEMALE', 'ILEUM_CONTROL_MALE', 'ILEUM_VEHICLE_FEMALE', 'ILEUM_VEHICLE_MALE', 'ILEUM_PROBIOTIC_FEMALE', 'ILEUM_PROBIOTIC_MALE', 'ILEUM_ANTIBIOTIC_FEMALE', 'ILEUM_ANTIBIOTIC_MALE', 'CECUM_CONTROL_FEMALE', 'CECUM_CONTROL_MALE', 'CECUM_VEHICLE_FEMALE', 'CECUM_VEHICLE_MALE', 'CECUM_PROBIOTIC_FEMALE', 'CECUM_PROBIOTIC_MALE', 'CECUM_ANTIBIOTIC_FEMALE', 'CECUM_ANTIBIOTIC_MALE' ) sample.item.color <- c( "#A9CCE3","#A9CCE3", "#2980B9","#2980B9", "#1B4F72","#1B4F72", "#001B3A","#001B3A", "#D2B4DE","#D2B4DE", "#8E44AD","#8E44AD", "#4A235A","#4A235A", "#1C0028","#1C0028", "#ABEBC6","#ABEBC6", "#2ECC71","#2ECC71", "#186A3B","#186A3B", "#144414","#144414", "#EDBB99","#EDBB99", "#D35400","#D35400", "#6E2C00","#6E2C00", "#431E03","#431E03" ) sample.item.linetype <- c( 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid', 'dashed','solid' ) ######## ### FULL ######## rareres <- ggrarecurve.get_rarecurve(project.ps, step=100, file=paste0(caminho,'Rarecurve_DADA2','.svg') ) svg(filename=paste0(caminho,'Rarecurve_SEX_DADA2','.svg') ) prare.sex <- ggrarecurve.rarecurve(obj=rareres, factorNames="SEX", indexNames=c("Observe"), shadow=TRUE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(sex.item.color, sex.item))+ scale_fill_manual(values=setNames(sex.item.color, sex.item))+ scale_linetype_manual(values=setNames(sex.item.linetype, sex.item))+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.sex graphics.off() svg(filename=paste0(caminho,'Rarecurve_TREATMENT_DADA2','.svg') ) prare.treatment <- ggrarecurve.rarecurve(obj=rareres, factorNames="TREATMENT", indexNames=c("Observe"), shadow=TRUE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(treatment.item.color, treatment.item), breaks=treatment.item )+ scale_fill_manual(values=setNames(treatment.item.color, treatment.item), breaks=treatment.item )+ scale_linetype_manual(values=setNames(treatment.item.linetype, treatment.item), breaks=treatment.item )+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.treatment graphics.off() svg(filename=paste0(caminho,'Rarecurve_SEGMENT_DADA2','.svg') ) prare.segment <- ggrarecurve.rarecurve(obj=rareres, factorNames="SEGMENT", indexNames=c("Observe"), shadow=TRUE, se=TRUE, linesize=1 ) + scale_color_manual(values=setNames(segment.item.color, segment.item), breaks=segment.item )+ scale_fill_manual(values=setNames(segment.item.color, segment.item), breaks=segment.item )+ scale_linetype_manual(values=setNames(segment.item.linetype, segment.item), breaks=segment.item )+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.segment graphics.off() svg(filename=paste0(caminho,'Rarecurve_SAMPLE_DADA2','.svg') ) prare.sample <- ggrarecurve.rarecurve(obj=rareres, factorNames=c("SAMPLE"), indexNames=c("Observe"), shadow=TRUE, se=TRUE, linesize=1, factorColor="SAMPLE", factorLine="SAMPLE" ) + scale_color_manual(values=setNames(sample.item.color, sample.item), breaks=sample.item )+ scale_fill_manual(values=setNames(sample.item.color, sample.item), breaks=sample.item )+ scale_linetype_manual(values=setNames(sample.item.linetype, sample.item), breaks=sample.item )+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.sample graphics.off() ######## ### RAREFIED ######## rareres.rarefied <- ggrarecurve.get_rarecurve(project.final.ps, step=100, file=paste0(caminho,'Rarecurve_DADA2_rarefied','.svg') ) svg(filename=paste0(caminho,'Rarecurve_SEX_DADA2_rarefied','.svg') ) prare.sex.rarefied <- ggrarecurve.rarecurve(obj=rareres.rarefied, factorNames="SEX", indexNames=c("Observe"), shadow=TRUE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(sex.item.color, sex.item))+ scale_fill_manual(values=setNames(sex.item.color, sex.item))+ scale_linetype_manual(values=setNames(sex.item.linetype, sex.item))+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.sex.rarefied graphics.off() svg(filename=paste0(caminho,'Rarecurve_TREATMENT_DADA2_rarefied','.svg') ) prare.treatment.rarefied <- ggrarecurve.rarecurve(obj=rareres.rarefied, factorNames="TREATMENT", indexNames=c("Observe"), shadow=TRUE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(treatment.item.color, treatment.item), breaks=treatment.item )+ scale_fill_manual(values=setNames(treatment.item.color, treatment.item), breaks=treatment.item )+ scale_linetype_manual(values=setNames(treatment.item.linetype, treatment.item), breaks=treatment.item )+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.treatment.rarefied graphics.off() svg(filename=paste0(caminho,'Rarecurve_SEGMENT_DADA2_rarefied','.svg') ) prare.segment.rarefied <- ggrarecurve.rarecurve(obj=rareres.rarefied, factorNames="SEGMENT", indexNames=c("Observe"), shadow=TRUE, se=TRUE, linesize=1 ) + scale_color_manual(values=setNames(segment.item.color, segment.item), breaks=segment.item )+ scale_fill_manual(values=setNames(segment.item.color, segment.item), breaks=segment.item )+ scale_linetype_manual(values=setNames(segment.item.linetype, segment.item), breaks=segment.item )+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.segment.rarefied graphics.off() svg(filename=paste0(caminho,'Rarecurve_SAMPLE_DADA2_rarefied','.svg') ) prare.sample.rarefied <- ggrarecurve.rarecurve(obj=rareres.rarefied, factorNames=c("SAMPLE"), indexNames=c("Observe"), shadow=TRUE, se=TRUE, linesize=1, factorColor="SAMPLE", factorLine="SAMPLE" ) + scale_color_manual(values=setNames(sample.item.color, sample.item), breaks=sample.item )+ scale_fill_manual(values=setNames(sample.item.color, sample.item), breaks=sample.item )+ scale_linetype_manual(values=setNames(sample.item.linetype, sample.item), breaks=sample.item )+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="grey"), strip.text.x = element_text(face="bold")) prare.sample.rarefied graphics.off() ``` ```{r, label="basic-rarefaction-curves-shown"} #| eval: true #| echo: false #| warning: false #| lightbox: true #| layout-nrow: 1 #| layout-ncol: 2 #| label: fig-rarecurve #| fig-cap: Curvas de Rarefação #| fig-subcap: #| - "Rarefaction curve (original counts)" #| - "Rarefaction curve (rarefied counts)" knitr::include_graphics(paste0(caminho,'Rarecurve_DADA2','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_DADA2_rarefied','.svg')) ``` ```{r, label="rarefaction-curves-shown"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 4 #| layout-nrow: 2 #| label: fig-rarecurve-groups #| fig-cap: Rarefaction analysis by groups #| fig-subcap: #| - "Rarefaction curve (original counts) by SEX" #| - "Rarefaction curve (original counts) by TREATMENT" #| - "Rarefaction curve (original counts) by SEGMENT" #| - "Rarefaction curve (original counts) by SAMPLE (SEGMENT+TREATMENT+SEX)" #| - "Rarefaction curve (rarefied counts) by SEX" #| - "Rarefaction curve (rarefied counts) by TREATMENT" #| - "Rarefaction curve (rarefied counts) by SEGMENT" #| - "Rarefaction curve (rarefied counts) by SAMPLE (SEGMENT+TREATMENT+SEX)" #| eval: true knitr::include_graphics(paste0(caminho,'Rarecurve_SEX_DADA2','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_TREATMENT_DADA2','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_SEGMENT_DADA2','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_SAMPLE_DADA2','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_SEX_DADA2_rarefied','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_TREATMENT_DADA2_rarefied','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_SEGMENT_DADA2_rarefied','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_SAMPLE_DADA2_rarefied','.svg')) ``` A seguir temos um resultado do quanto a matriz de ASVs após a rarefação é esparsa utilizando a função *summarize_phyloseq()* que sumariza os resultados do objeto phyloseq, desta vez após a rarefação (`r spa.rarefied`). #### Shared OTUs A seguir, comparamos as ASVs compartilhadas entre as amostras considerando os dados completos e posteriormente os dados após o procedimento de rarefação descrito na seção anterior. ```{r, label="shared-otus-analysis"} #| echo: false #| warnings: false #| include: false #| eval: false # Create a Empty DataFrame otu.perc = as.data.table(data.frame(matrix(nrow = dim(sample_data(project.ps))[1], ncol = dim(sample_data(project.ps))[1]) )) rownames(otu.perc) <- rownames(sample_data(project.ps)) colnames(otu.perc) <- rownames(sample_data(project.ps)) for (col in colnames(otu.perc)) { otu.perc[[col]] <- as.numeric(otu.perc[[col]]) } r <- setSharedPerc(otu.perc, project.ps ) otu.perc <- as.data.frame(x=otu.perc, row.names=rownames(otu.perc)) file.remove(paste0(caminho,'Shared-ASVs.xlsx')) wb <- openxlsx::createWorkbook() openxlsx::addWorksheet(wb, "Shared ASVs (original)") openxlsx::writeData(wb, "Shared ASVs (original)", otu.perc, startRow = 1, startCol = 1) otu.rare.perc = as.data.table(data.frame(matrix(nrow = dim(sample_data(project.final.ps))[1], ncol = dim(sample_data(project.final.ps))[1]) )) rownames(otu.rare.perc) <- rownames(sample_data(project.final.ps)) colnames(otu.rare.perc) <- rownames(sample_data(project.final.ps)) for (col in colnames(otu.rare.perc)) { otu.rare.perc[[col]] <- as.numeric(otu.rare.perc[[col]]) } r <- setSharedPerc(otu.rare.perc, project.final.ps) otu.rare.perc <- as.data.frame(x=otu.rare.perc, row.names=rownames(otu.rare.perc)) addWorksheet(wb, "Shared ASVs (rarefied)") writeData(wb, "Shared ASVs (rarefied)", otu.rare.perc, startRow = 1, startCol = 1) saveWorkbook(wb, file=paste0(caminho,'Shared-ASVs.xlsx') , overwrite = TRUE) otu.perc.sample.annotation = data.frame(SEGMENT=as.factor(sample_data(project.ps)[['SEGMENT']]), TREATMENT=as.factor(sample_data(project.ps)[['TREATMENT']]), SEX=as.factor(sample_data(project.ps)[['SEX']]) ) svg(filename=paste0(caminho,'aHeatmap_shared_otus_original.svg'),width = 20,height=30 ) aheatmap(data.matrix(otu.perc), color = 1L, Rowv=NA, Colv=NA, annCol=otu.perc.sample.annotation, annRow=otu.perc.sample.annotation, annColors = list(SEGMENT=setNames(segment.item.color, as.character(segment.item))[levels(otu.perc.sample.annotation$SEGMENT)], TREATMENT=setNames(treatment.item.color, as.character(treatment.item))[levels(otu.perc.sample.annotation$TREATMENT)], SEX=setNames(sex.item.color, as.character(sex.item))[levels(otu.perc.sample.annotation$SEX)] ) ) graphics.off() svg(filename=paste0(caminho,'aHeatmap_shared_otus_rarefied.svg'),width = 20,height=30 ) aheatmap(data.matrix(otu.rare.perc), color = 1L, Rowv=NA, Colv=NA, annCol=otu.perc.sample.annotation, annRow=otu.perc.sample.annotation, annColors = list(SEGMENT=setNames(segment.item.color, as.character(segment.item))[levels(otu.perc.sample.annotation$SEGMENT)], TREATMENT=setNames(treatment.item.color, as.character(treatment.item))[levels(otu.perc.sample.annotation$TREATMENT)], SEX=setNames(sex.item.color, as.character(sex.item))[levels(otu.perc.sample.annotation$SEX)] ) ) graphics.off() ``` ```{r, label="aheatmap-shared-otus"} #| echo: false #| eval: true #| warning: false #| lightbox: true #| layout-ncol: 2 #| layout-nrow: 1 #| label: fig-shared-otus-groups #| fig-cap: Rarefaction analysis by groups #| fig-subcap: #| - "Percentage of shared ASVs (original dataset) by GROUP" #| - "Percentage of shared ASvs (rarefied dataset) by GROUP" knitr::include_graphics(paste0(caminho,'aHeatmap_shared_otus_original','.svg')) knitr::include_graphics(paste0(caminho,'aHeatmap_shared_otus_rarefied','.svg')) ``` #### Aglomeração de ASVs da mesma taxonomia ```{r, label="taxonomy-agglomeration"} #| echo: false #| warning: false #| eval: false #| include: false project.rarefied.genus.ps = tax_glom(project.final.ps, taxrank="Genus", NArm=FALSE) ``` Após a utilização da função *tax_glom()* selecionando o nível de gênero, na tabela de ASVs temos `r dim(otu_table(project.rarefied.genus.ps))[1]` registros. #### Análise de intersecção As análises de intersecção a seguir estão organizadas por grupos de amostras (SEGMENTO), portanto, os números relacionados ao grau zero de intersecção são as que não entraram na comparação, estão nos demais grupos. ```{r, label="upset-plots"} #| echo: false #| warning: false #| include: false #| eval: false metadata$TREATMENT <- as.factor(as.character(metadata$TREATMENT)) metadata$SEGMENT <- as.factor(as.character(metadata$SEGMENT)) metadata$SAMPLE <- as.factor(as.character(metadata$SAMPLE)) metadata$SEX <- as.factor(as.character(metadata$SEX)) for (seg in unique(metadata$SEGMENT)) { metadata.seg <- subset(metadata, SEGMENT==seg) upsetda <- get_upset(obj=as.data.frame(t(as.matrix(otu_table(project.final.ps)[,rownames(metadata.seg)]))), sampleda=metadata.seg, factorNames='TREATMENT') m1 = make_comb_mat(upsetda, mode = "intersect",top_n_sets = 20) upset.asvs<-UpSet(m1, right_annotation = upset_right_annotation(m1, add_numbers = TRUE), top_annotation = HeatmapAnnotation( degree = as.character(comb_degree(m1)), col = list(degree = setNames(rev(brewer.pal(length(unique(comb_degree(m1))),"YlGn")), as.character(sort(unique(comb_degree(m1)),decreasing=TRUE)) ) ), "Intersection\nsize" = anno_barplot(comb_size(m1), border = FALSE, add_numbers = TRUE, gp = gpar(fill = "black"), height = unit(2, "cm") ), annotation_name_side = "left", annotation_name_rot = 0) ) svg(paste0(caminho,'UpSet_asvs_',seg,'.svg') ) plot(upset.asvs) graphics.off() upsetda.genus <- get_upset(obj=as.data.frame(t(as.matrix(otu_table(project.rarefied.genus.ps)[,rownames(metadata.seg)]))), sampleda=metadata.seg, factorNames='TREATMENT') m2 = make_comb_mat(upsetda.genus, mode = "intersect",top_n_sets = 20) upset.genus <-UpSet(m2, right_annotation = upset_right_annotation(m2, add_numbers = TRUE), top_annotation = HeatmapAnnotation( degree = as.character(comb_degree(m2)), col = list(degree = setNames(rev(brewer.pal(length(unique(comb_degree(m2))),"YlGn")), as.character(sort(unique(comb_degree(m2)),decreasing=TRUE)) ) ), "Intersection\nsize" = anno_barplot(comb_size(m2), border = FALSE, add_numbers = TRUE, gp = gpar(fill = "black"), height = unit(2, "cm") ), annotation_name_side = "left", annotation_name_rot = 0) ) svg(paste0(caminho,'UpSet_genus_',seg,'.svg') ) plot(upset.genus) graphics.off() } ``` ```{r, label="upset-plots-by_group-shown"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 2 #| layout-nrow: 4 #| label: fig-upset #| fig-cap: Intersection analyses (UpSet plots) among TREATMENTS by SEGMENT #| fig-subcap: #| - "ASV level (DUODENUM)" #| - "Genus level (DUODENUM)" #| - "ASV level (JEJUNUM)" #| - "Genus level (JEJUNUM)" #| - "ASV level (ILEUM)" #| - "Genus level (ILEUM)" #| - "ASV level (CECUM)" #| - "Genus level (CECUM)" #| eval: true knitr::include_graphics(paste0(caminho,'UpSet_asvs_','DUODENUM','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_genus_','DUODENUM','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_asvs_','JEJUNUM','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_genus_','JEJUNUM','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_asvs_','ILEUM','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_genus_','ILEUM','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_asvs_','CECUM','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_genus_','CECUM','.svg')) ``` #### Composição taxonômica A seguir temos a composição dos mais abundantes filos e gêneros em cada amostras (@fig-compo). Os gráficos foram gerados a partir do pacote microeco v1.6.0. Para isso, o objeto phyloseq foi convertido para um objeto microeco por meio da função *phyloseq2meco()* do pacote file2meco v0.7.1. ```{r, label="taxonomy-composition"} #| echo: false #| warning: false #| include: false #| eval: false meco <- phyloseq2meco(project.final.ps) meco$sample_table$SEX <- c('M','F')[ match(meco$sample_table$SEX, c('MALE','FEMALE') ) ] #plot_bar(project.rarefied.genus.ps, fill = "Genus") # hcl.colors(10, palette = "Earth") colorsEarth = c("#A36B2B", "#B2874C", "#C2A36C", "#D2C08D", "#E4DCB0", "#D8DFC3", "#B7C7B8", "#87B4A3", "#509F9E", "#2686A0") colorsEarth = hcl.colors(10, palette = "Earth") # hcl.colors(10, palette = "Fall") colorsFall = c("#3C5941", "#61795B", "#899B76", "#B4BD93", "#E2E0B3", "#F2E2B3", "#E3C28D", "#D9A066", "#D07C42", "#C7522B") colorsFall = hcl.colors(5, palette = "Fall") compPhylum <- trans_abund$new(dataset = meco, taxrank = "Phylum", ntaxa = 5) phylumPlot <- compPhylum$plot_bar(others_color = "grey90", facet = c("SEGMENT","TREATMENT","SEX"), xtext_keep = TRUE, xtext_angle = 45, legend_text_italic = FALSE, color_values = colorsFall)+theme(axis.text.x=element_text(angle=90,vjust=0.5)) svg(paste0(caminho,'CompoPlot_','phylum','.svg'), width=20 ) phylumPlot graphics.off() compGenus <- trans_abund$new(dataset = meco, taxrank = "Genus", ntaxa = 10) genusPlot <- compGenus$plot_bar(others_color = "grey90", facet = c("SEGMENT","TREATMENT","SEX"), xtext_keep = TRUE, legend_text_italic = FALSE, color_values = colorsEarth)+theme(axis.text.x=element_text(angle=90,vjust=0.5)) svg(paste0(caminho,'CompoPlot_','genus','.svg'), width=20 ) genusPlot graphics.off() ``` ```{r, label="taxonomy-composition-shown"} #| echo: false #| warning: false #| lightbox: true #| layout-nrow: 2 #| layout-ncol: 1 #| label: fig-compo #| fig-cap: Compositional plot #| fig-subcap: #| - "Phylum level (Top 5)" #| - "Genus level (Top 10)" #| eval: true knitr::include_graphics(paste0(caminho,'CompoPlot_phylum','.svg')) knitr::include_graphics(paste0(caminho,'CompoPlot_genus','.svg')) ``` #### Alfa diversidade ```{r, label="alpha-diversity-indices"} #| echo: false #| warning: false #| eval: false #alpha.indices <- c("Observed", "Chao1", "ACE", "Simpson", "InvSimpson", "Shannon", "Fisher") #div.table <- estimate_richness(project.final.ps, split=TRUE, measures=alpha.indices) # https://mibwurrepo.github.io/Microbial-bioinformatics-introductory-course-Material-2018/alpha-diversities.html div.table <- alpha(otu_table(project.final.ps, taxa_are_rows=TRUE), index="all") div.table[['dominance_core_abundance']] <- core_abundance(otu_table(project.final.ps, taxa_are_rows=TRUE), detection = 0.1/100, prevalence = 45/100) df.pd <- pd(t(as.matrix(otu_table(project.final.ps, taxa_are_rows=TRUE))), phy_tree(project.final.ps),include.root=T) div.table[rownames(div.table),'PD'] <- df.pd[rownames(div.table),'PD'] openxlsx::write.xlsx(x=div.table, file=paste0(caminho,'diversity.xlsx'), sheetName = "Alpha Diversity Indices") ``` Os índices de alfa diversidade (`r str_escape( paste(gsub('(.*)','"\\1"',as.character( colnames(div.table) )),collapse=", "))`) foram obtidos com a função *alpha()* do pacote microbiome v.1.16.0. No resultado, há métricas para avaliar riqueza, diversidade, equitabilidade, dominância e raridade. Além disso, também foi calculada uma métrica de diversidade filogenética (Faith's Phylogenetic Diversity), calculada pela função *pd()* pacote *picante* v1.8.2. O resultado (@tbl-alphadiversity) foi incorporado também à planilha "Alpha Diversity Indices" do Excel (@sec-planilha). Para informações sobre os índices de diversidade, incluindo a forma como interpretar os resultados, consulte a seção "[Indices of diversity and eveness](https://www.davidzeleny.net/anadat-r/doku.php/en:div-ind)" da página "[Analysis of community ecology data in R](https://www.davidzeleny.net/anadat-r/doku.php/en:start)" de David Zelený". ##### Índices de Alfa diversidade - Riqueza (referem-se à quantidade de diferentes taxa): - observed; - chao1; - Diversidade: - inverse_simpson; - gini_simpson; - shannon; - fisher; - coverage (refere-se ao número de grupos necessários para terem uma determinada proporção no ecossistema, aqui foi utilizado o valor pré-determinado de 0,5, *i.e.* 50%); - PD: Faith's Phylogenetic Diversity (calcula a soma do comprimento das ramificações dentro de uma amostra); - Dominância (referem-se à abundância dos taxa mais abundantes): - dbp; - dmn; - absolute; - relative; - simpson; - core_abundance (refere-se à proporção relativa dos taxa essenciais "core", *i.e.* função *core_abundance()* com os parâmetros: detection = .1/100 e prevalence = 45/100); - gini; - Raridade e baixa abundância (quantificam a concentração de taxa raros ou pouco abundantes): - log_modulo_skewness; - low_abundance; - rare_abundance (complemento de core_abundance: 1-*core_abundance()*, considerando os parâmetros: detection = 0.2/100, prevalence = 20/100); - Equitabilidade: - camargo; - pielou; - simpson; - evar; - bulla; ```{r, label="alpha-diversity-indices-shown"} #| echo: false #| warning: false #| eval: true xfun::embed_file(paste0(caminho,'diversity.xlsx')) ``` ```{r, label="alpha-diversity-indices-table"} #| echo: false #| label: tbl-alphadiversity #| tbl-cap: Alpha diversity indices #| eval: true knitr::kable(div.table) ``` ##### Comparações de Alfa diversidade Também foram realizadas comparações entre as crianças considerando um índice de riqueza (observed) e um outro de diversidade (shannon). A @fig-plotalphadiv exibe as diferenças referente a esses índices nas amotras por segmento do intestino e tratamento (SEGMENT_TREATMENT). Considerando cada amostra um ambiente, podemos considerar que a diversidade relacionada ao conjunto de todas as amostras (POOLs) como sendo a diversidade gama (@fig-plotgammadiv), ou seja, são os mesmos índices de diversidade alfa porém considerando todo o conjunto de amostras consideradas no sequenciamento. ```{r, label="alpha_and_gamma-div-plot"} #| echo: false #| warning: false #| include: false #| eval: false metadata.div <- metadata metadata.div$SEGMENT_TREATMENT <- paste0(metadata.div$SEGMENT, '_', metadata.div$TREATMENT) metadata.div$SEQUENCING <- 1 for (i in colnames(div.table)) { metadata.div[[i]] <- div.table[[i]] } metadata.div$SEGMENT_TREATMENT <- as.factor(as.character(metadata.div$SEGMENT_TREATMENT)) # create a list of pairwise comaprisons cmp <- levels(metadata.div$SEGMENT_TREATMENT) # get the variables # make a pairwise list that we want to compare. cmp.pairs <- combn(seq_along(cmp), 2, simplify = FALSE, FUN = function(i)cmp[i]) palette <- distinctColorPalette(length(cmp)) project.cols <- setNames(palette,cmp) for (i in colnames(div.table)) { p1 <- ggviolin(metadata.div, x = "TREATMENT", y = i, add = "boxplot", fill = "TREATMENT") + scale_fill_manual(values=setNames(treatment.item.color, treatment.item), breaks=treatment.item) p1 <- p1 + geom_pwc(method="wilcox_test", label = "{p.format} {p.signif}", hide.ns = 'p') + facet_wrap(~SEGMENT)+ scale_x_discrete(guide = guide_axis(angle = 45)) ggsave(paste0(caminho,'AlphaDivPlot_',i,'.svg'), p1, width=20, unit='cm' ) p2 <- ggviolin(metadata.div, x = "SEQUENCING", y = i, add = "boxplot", fill = "gray") ggsave(paste0(caminho,'GammaDivPlot_',i,'.svg'), p2, width=20, unit='cm' ) } ``` ```{r, label="alpha-div-plot-shown"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 2 #| label: fig-plotalphadiv #| fig-cap: Alpha diversity indices by segment and treatment #| fig-subcap: #| - "Richness index (observed)" #| - "Diversity index (shannon)" #| eval: true knitr::include_graphics(paste0(caminho,'AlphaDivPlot_observed','.svg')) knitr::include_graphics(paste0(caminho,'AlphaDivPlot_diversity_shannon','.svg')) ``` ```{r, label="gamma-div-plot-shown"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 2 #| label: fig-plotgammadiv #| fig-cap: Gamma diversity indices #| fig-subcap: #| - "Richness index (observed)" #| - "Diversity index (shannon)" #| eval: true knitr::include_graphics(paste0(caminho,'GammaDivPlot_observed','.svg')) knitr::include_graphics(paste0(caminho,'GammaDivPlot_diversity_shannon','.svg')) ``` Um outro método para avaliação da alfa diversidade é a baseada nos números de Hill, portanto, realizamos uma análise com o pacote *hilldiv* v1.5.1 do R (@fig-hillprofile). ```{r, label="hill-profile"} #| echo: false #| warning: false #| include: false #| eval: false metadata$SEGMENT_TREATMENT <- paste0(metadata$SEGMENT, '_', metadata$TREATMENT) ultrametric_tree <- force.ultrametric( phy_tree(project.final.ps) ) hierarchy <- data.frame(Sample=rownames(sample_data(project.final.ps)), Group=metadata[rownames(sample_data(project.final.ps)),'SAMPLE']) svg(paste0(caminho,'Hill_profile_alpha_SAMPLE','.svg') ) div_profile_plot(div_profile(count=otu_table(project.final.ps), hierarchy=hierarchy ),colour=setNames(sample.item.color, sample.item)) graphics.off() hierarchy <- data.frame(Sample=rownames(sample_data(project.final.ps)), Group=metadata[rownames(sample_data(project.final.ps)),'SEGMENT_TREATMENT']) svg(paste0(caminho,'Hill_profile_alpha_SEGMENT_TREATMENT','.svg') ) ### div_profile_plot(div_profile(count=otu_table(project.final.ps), ### tree=ultrametric_tree ### )) segment_treatment.item <- unique(gsub('_(MALE|FEMALE)$','',sample.item)) segment_treatment.item.color <- unique(sample.item.color) div_profile_plot(div_profile(count=otu_table(project.final.ps), hierarchy=hierarchy ), colour=setNames(segment_treatment.item.color, segment_treatment.item), log=FALSE) graphics.off() svg(paste0(caminho,'Hill_profile_alpha_log_SEGMENT_TREATMENT','.svg') ) div_profile_plot(div_profile(count=otu_table(project.final.ps), hierarchy=hierarchy ), colour=setNames(segment_treatment.item.color, segment_treatment.item), log=TRUE) graphics.off() svg(filename=paste0(caminho,'Hill_profile_gamma','.svg')) div_profile_plot(div_profile(count=otu_table(project.final.ps), level="gamma" )) graphics.off() ``` ```{r, label="hill-profile-shown"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 4 #| label: fig-hillprofile #| fig-cap: Hill diversity plots #| fig-subcap: #| - "Alpha diversity profiles by SAMPLE" #| - "Alpha diversity profiles by SEGMENT & TREATMENT" #| - "Alpha diversity profiles by SEGMENT & TREATMENT (Log-transformed)" #| - "Gamma diversity profile" #| eval: true knitr::include_graphics(paste0(caminho,'Hill_profile_alpha_SAMPLE','.svg')) knitr::include_graphics(paste0(caminho,'Hill_profile_alpha_SEGMENT_TREATMENT','.svg')) knitr::include_graphics(paste0(caminho,'Hill_profile_alpha_log_SEGMENT_TREATMENT','.svg')) knitr::include_graphics(paste0(caminho,'Hill_profile_gamma','.svg')) ``` #### Beta diversidade Também foram realizadas análises de Beta diversidade comparando as amostras usando Análises de Coordenadas Principais (PCoA) com distâncias categóricas de Jaccard (presença/ausência) e Bray-curtis (quantitativa), e as distâncias UniFrac que levam em conta as relações filogenéticas. As distâncias UniFrac consideradas foram a ponderada, em que as diferenças entre as abundâncias dos *taxa* são consideradas, e a não-ponderada, em que não se considera essas diferenças de abundâncias. Veja apresentação sobre distâncias neste [link](https://wiki.qcbs.ca/_media/lecture6a_distances_ordination.pdf) e métodos de ordenação neste [link](https://xratio.medium.com/pca-pcoa-mds-fa-ca-demystify-dimensionality-reduction-techniques-2-8a65d79c2e8b). A diversidade beta observada pode ser avaliada por meio das figuras a seguir (@fig-betadiv-sample e @fig-betadiv-segment_treatment). ```{r, label="beta-div-plot-1"} #| echo: false #| warning: false #| include: false #| eval: false meco$sample_table[rownames(meco$sample_table),'SEGMENT_TREATMENT'] <- metadata[rownames(meco$sample_table),'SEGMENT_TREATMENT'] meco$cal_betadiv(unifrac = TRUE) t1 <- trans_beta$new(dataset = meco, group = "SEGMENT_TREATMENT", measure = "bray") t1$ordination_method <- "PCoA" t1$cal_ordination() ggsave(paste0(caminho,'PCoA_bray_SEGMENT_TREATMENT','.svg'), t1$plot_ordination(plot_color="SEGMENT_TREATMENT", color_values = setNames(segment_treatment.item.color, segment_treatment.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) t2 <- trans_beta$new(dataset = meco, group = "SEGMENT_TREATMENT", measure = "jaccard") t2$ordination_method <- "PCoA" t2$cal_ordination() ggsave(paste0(caminho,'PCoA_jaccard_SEGMENT_TREATMENT','.svg'), t2$plot_ordination(plot_color="SEGMENT_TREATMENT", color_values = setNames(segment_treatment.item.color, segment_treatment.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) t3 <- trans_beta$new(dataset = meco, group = "SEGMENT_TREATMENT", measure = "wei_unifrac") t3$ordination_method <- "PCoA" t3$cal_ordination() ggsave(paste0(caminho,'PCoA_wei_unifrac_SEGMENT_TREATMENT','.svg'), t3$plot_ordination(plot_color="SEGMENT_TREATMENT", color_values = setNames(segment_treatment.item.color, segment_treatment.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) t4 <- trans_beta$new(dataset = meco, group = "SEGMENT_TREATMENT", measure = "unwei_unifrac") t4$ordination_method <- "PCoA" t4$cal_ordination() ggsave(paste0(caminho,'PCoA_unwei_unifrac_SEGMENT_TREATMENT','.svg'), t4$plot_ordination(plot_color="SEGMENT_TREATMENT", color_values = setNames(segment_treatment.item.color, segment_treatment.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) ``` ```{r, label="beta-div-plot-2"} #| echo: false #| warning: false #| include: false #| eval: false t1 <- trans_beta$new(dataset = meco, group = "SAMPLE", measure = "bray") t1$ordination_method <- "PCoA" t1$cal_ordination() ggsave(paste0(caminho,'PCoA_bray_SAMPLE','.svg'), t1$plot_ordination(plot_color="SEGMENT_TREATMENT", color_values = setNames(sample.item.color, sample.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) t2 <- trans_beta$new(dataset = meco, group = "SAMPLE", measure = "jaccard") t2$ordination_method <- "PCoA" t2$cal_ordination() ggsave(paste0(caminho,'PCoA_jaccard_SAMPLE','.svg'), t2$plot_ordination(plot_color="SAMPLE", color_values = setNames(sample.item.color, sample.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) t3 <- trans_beta$new(dataset = meco, group = "SAMPLE", measure = "wei_unifrac") t3$ordination_method <- "PCoA" t3$cal_ordination() ggsave(paste0(caminho,'PCoA_wei_unifrac_SAMPLE','.svg'), t3$plot_ordination(plot_color="SAMPLE", color_values = setNames(sample.item.color, sample.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) t4 <- trans_beta$new(dataset = meco, group = "SAMPLE", measure = "unwei_unifrac") t4$ordination_method <- "PCoA" t4$cal_ordination() ggsave(paste0(caminho,'PCoA_unwei_unifrac_SAMPLE','.svg'), t4$plot_ordination(plot_color="SAMPLE", color_values = setNames(sample.item.color, sample.item), plot_type = c("point", "ellipse")), width=20, unit='cm' ) ``` ```{r, label="beta-div-plot-shown-2"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 2 #| layout-nrow: 2 #| label: fig-betadiv-sample #| fig-cap: Beta diversity analyses using PCoA ordination (SAMPLE) #| fig-subcap: #| - "Bray-curtis distance" #| - "Jaccard distance" #| - "Weighted UniFrac distance" #| - "Unweighted UniFrac distance" #| eval: true knitr::include_graphics(paste0(caminho,'PCoA_bray_SAMPLE','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_jaccard_SAMPLE','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_wei_unifrac_SAMPLE','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_unwei_unifrac_SAMPLE','.svg')) ``` ```{r, label="beta-div-plot-shown-1"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 2 #| layout-nrow: 2 #| label: fig-betadiv-segment_treatment #| fig-cap: Beta diversity analyses using PCoA ordination (SEGMENT_TREATMENT) #| fig-subcap: #| - "Bray-curtis distance" #| - "Jaccard distance" #| - "Weighted UniFrac distance" #| - "Unweighted UniFrac distance" #| eval: true knitr::include_graphics(paste0(caminho,'PCoA_bray_SEGMENT_TREATMENT','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_jaccard_SEGMENT_TREATMENT','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_wei_unifrac_SEGMENT_TREATMENT','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_unwei_unifrac_SEGMENT_TREATMENT','.svg')) ``` ### Análise de abundância diferencial Aqui podemos fazer uma comparação da abundância dos *taxa* (gêneros) entre as microbiotas nos segmentos do intestino sob um determinado tratamento, para os dois pools de amostras de aves, considerando machos e fêmeas. A avaliação da abundância diferencial foi realizada com o pacote DESeq2 v1.34.0, a partir da matrix após a filtragem dos *taxa* com valores de abundância baixos de acordo com a abundância mínima de 10 em ao menos 1 amostra. A normalização dos dados foi realizada seguida de transformação logarítmica. A abundância diferencial foi avaliada aplicando modelo linear generalizado (*Generalized Linear Model* - GLM) considerando um ajuste a uma distribuição de probabilidade binomial negativa (*Negative Binomial GLM*). As significâncias dos coeficientes do GLM (gêneros) foi avaliada com o teste de Wald (*Generalized Linear Model Likelihood Ratio Test*). Os dados tiveram seus níveis descritivos de significância estatística (*p-values*) ajustados (*adjusted p-values*) pela taxa de falsas descobertas (*False Discovery Rate*, ou FDR). Os dados de abundância das amostras exibidos e suas médias foram previamente transformados utlizando a transformação Log2 (função *normTransform()* do pacote DESeq2 com parâmetro *f=log2* e *pc=1* para indicar pseudocontagem de 1). ```{r, label="differential-abundance"} #| echo: false #| warning: false #| include: false #| eval: false counts.deg <- as(phyloseq::otu_table(project.rarefied.genus.ps, taxa_are_rows=TRUE), "matrix") lineage <- apply((as(phyloseq::tax_table(project.rarefied.genus.ps),"matrix")), 1, function(x) { return(paste(x[c('Kingdom','Phylum','Class','Order','Family','Genus')],collapse=";")) } ) rownames(counts.deg) <- as.character(lineage[rownames(counts.deg)]) metadata.deg <- metadata for (col in c('PREFIX','SAMPLE','POOL','TREATMENT','SEX','SEGMENT','SEGMENT_TREATMENT')) { metadata.deg[[col]] <- as.factor(as.character(metadata.deg[[col]])) } metadata.deg$TREATMENT <- relevel(metadata.deg$TREATMENT, ref="CONTROL") dds <- DESeqDataSetFromMatrix(countData = counts.deg, colData = DataFrame( metadata.deg[colnames(counts.deg),] ), design = ~0+SEGMENT:TREATMENT, tidy=FALSE) #model.matrix(~0+SEGMENT+TREATMENT+SEX+SEGMENT:TREATMENT, data=DataFrame( metadata.deg[colnames(counts.deg),] ) ) smallestGroupSize <- 1 smallestLibSize <- 10 smallestCount <- 10 keep_rows <- rowSums(counts(dds) >= smallestCount) >= smallestGroupSize keep_cols <- colSums(counts(dds)) >= smallestLibSize excLibs <- names( colSums(counts(dds))[colSums(counts(dds)) <= smallestLibSize] ) #names(keep_rows[keep_rows]) dds <- dds[keep_rows, keep_cols] #dds <- estimateSizeFactors(dds) #dds <- estimateDispersions(dds) #dds <- nbinomWaldTest(dds) dds <- DESeq(dds) ddsClean <- dds[which(mcols(dds)$betaConv),] ntd <- normTransform(ddsClean) #vsd <- vst(dds, blind=FALSE) #vsd <- varianceStabilizingTransformation(counts(dds)) #rld <- rlog(dds, blind=FALSE) #meanSdPlot(assay(ntd)) #meanSdPlot(vsd) #meanSdPlot(assay(rld)) counts.ntd <- as.data.frame(assay(ntd)) prefix.treatment <- list('ANTIBIOTIC'='A', 'PROBIOTIC'='P', 'VEHICLE'='V', 'CONTROL'='C') letras <- list() file.remove(paste0(caminho,'Differential_Abundance.xlsx')) wb <- createWorkbook() for (s in levels(ddsClean$SEGMENT) ) { smpst <- rownames(subset(metadata.deg[names(keep_cols[keep_cols]),], SEGMENT==s )) counts.ntd[[paste0(s,'.','sum')]] <- rowSums(counts.ntd[,smpst]) counts.ntd[[paste0(s,'.','mean')]] <- rowSums(counts.ntd[,smpst])/length(smpst) for (t in levels(ddsClean$TREATMENT)) { smpst <- rownames(subset(metadata.deg, SEGMENT==s & TREATMENT==t)) if (length(smpst)>1) { counts.ntd[[paste0(s,'.','sum')]] <- rowSums(counts.ntd[,smpst]) counts.ntd[[paste0(s,'.',t,'.','mean')]] <- rowSums(counts.ntd[,smpst])/length(smpst) } else { counts.ntd[[paste0(s,'.','sum')]] <- counts.ntd[,smpst] counts.ntd[[paste0(s,'.',t,'.','mean')]] <- counts.ntd[,smpst]/length(smpst) } } all.comb.s <- apply(t(combn(paste0('SEGMENT',s,'.','TREATMENT',levels( metadata.deg$TREATMENT) ),2)), 1, function(x) { paste( x ,collapse="-") }) res <- list() compara <- data.frame(row.names=rownames( counts.ntd )) for (c in 1:length(all.comb.s)) { ctrst <- unlist(strsplit(all.comb.s[c],split='-')) new.comb.s <- gsub(paste0('SEGMENT',s,'.TREATMENT'),'',all.comb.s[c]) res[[ new.comb.s ]] <- results(ddsClean, list(ctrst[1],ctrst[2])) dalrt.df <- as.data.frame( res[[ new.comb.s ]] ) t1 <- gsub('^.*TREATMENT','',ctrst[1]) t2 <- gsub('^.*TREATMENT','',ctrst[2]) smps1 <- rownames(subset(metadata.deg, SEGMENT==s & TREATMENT==t1)) smps2 <- rownames(subset(metadata.deg, SEGMENT==s & TREATMENT==t2)) da.res <- (merge(x=dalrt.df, y=counts.ntd[,c(smps1, smps2, paste(s,t1,'mean',sep='.'), smps2, paste(s,t2,'mean',sep='.') )], by='row.names'))[,c('Row.names', smps1, paste(s,t1,'mean',sep='.'), smps2, paste(s,t2,'mean',sep='.'), 'log2FoldChange', 'pvalue', 'padj' )] rownames(da.res) <- da.res$Row.names wb.name <- paste0("DA_",s,'_',paste(c(as.character(prefix.treatment[t1]),as.character(prefix.treatment[t2])),collapse='-')) addWorksheet(wb, wb.name) writeData(wb, wb.name, da.res, startRow = 1, startCol = 1) #openxlsx::write.xlsx(x=da.res, file=paste0(caminho,'Differential_Abundance.xlsx'), sheetName = paste0("DA_",s,'_',paste(c(as.character(prefix.treatment[t1]),as.character(prefix.treatment[t2])),collapse='-')),overwrite=TRUE, append=TRUE) for (g in rownames(compara)) { compara[g,new.comb.s] <- ifelse(is.na(da.res[g,'padj']), FALSE, ifelse(da.res[g,'padj']<=0.05, TRUE, FALSE)) } } letras[[s]] <- data.frame(row.names=rownames(counts.ntd)) for (g in rownames(letras[[s]])) { difcompara <- as.logical(compara[g, colnames(compara) ]) names(difcompara) <- colnames(compara) difcomparaL <- multcompLetters(difcompara)$Letters letras[[s]][g,levels(ddsClean$TREATMENT)] <- difcomparaL[levels(ddsClean$TREATMENT)] } } counts.ntd[['sum']] <- rowSums(counts.ntd) counts.deg <- as.data.frame(counts.deg) counts.deg[['sum']] <- rowSums(counts.deg) # #(counts.deg[order(counts.deg$sum,decreasing=TRUE),])[1:30,levels(group)] saveWorkbook(wb, file=paste0(caminho,'Differential_Abundance.xlsx') , overwrite = TRUE) ``` Essa comparação foi realizada utilizando como entrada a aglomeração dos ASVs em gêneros e foram filtrados `r length(keep_rows[keep_rows==FALSE])` gêneros, sem que os grupos de amostras tivessem a soma das abundâncias mínima de `r smallestCount`. Assim como bibliotecas sem um mínimo de soma de abundâncias (`r smallestLibSize`), também foram removidas (`r paste(excLibs, collapse=".")`). Um Heatmap foi montado com os gêneros que tiveram alguma diferença de abundância significativa (**adj-p-value** \< 0.1) a partir da comparação de abundância diferencial entre os grupos em cada tempo (@fig-heamapda). Os dendrogramas foram construídos a partir do algoritmo de agrupamento hierárquico com distância *euclideana* e método de ligação *ward.D2*. Os valores para o gráfico foram escalados para z-score. Os valores relacionados ao cálculo de diferença de abundância foram inseridos em planilhas do arquivo Excel a seguir: ```{r, label="differential-abundance-shown"} #| echo: false #| warning: false #| eval: true xfun::embed_file(paste0(caminho,'Differential_Abundance.xlsx')) ``` ```{r, label="heatmap"} #| echo: false #| warning: false #| include: false #| eval: false # D | Duodenum | # J | Jejunum | # I | Ileum | # C | Cecum | for (s in segment.item) { cols <- paste0(s,'.',levels(ddsClean$TREATMENT),'.mean') any.signif <- apply(letras[[s]], 1, function(x) { length(unique(x))>1 } ) sum.ordered <- order(counts.ntd[[paste0(s,'.sum')]], decreasing=TRUE)[1:30] toplot <- t(scale(t(as.matrix( counts.ntd[any.signif,cols] )))) letras_toplot <- letras[[s]][rownames(toplot),] # letras_toplot <- letras[[s]][rownames(toplot),colnames(toplot)] col_fun <- colorRamp2(c(-4, 0, 4), c("green", "white", "red")) ht<-Heatmap(toplot, border=TRUE, col=col_fun, border_gp = gpar(col = "white"), rect_gp = gpar(col = "white", lwd = 2), show_heatmap_legend = TRUE, cluster_rows = TRUE, cluster_columns = TRUE, heatmap_legend_param = list(title_position = "leftcenter-rot", title = 'Z-score (abundance)'), cell_fun = function(j, i, x, y, w, h, col) { # since grid.text can also be vectorized grid.text(letras_toplot[i, j], x, y, gp = gpar(fontsize = 16)) } ) svg(filename=paste0(caminho,'Heatmap_DA_DESeq2_',s,'.svg'),width = 20,height=25 ) #bottom, left, top and right sides draw(ht, merge_legend = TRUE, heatmap_legend_side = "left", padding = unit(c(2, 2, 2, 20), 'cm')) graphics.off() } ``` ```{r, label="heatmap-shown"} #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 2 #| layout-nrow: 2 #| label: fig-heamapda #| fig-cap: Heatmap of top 30 most abundant genus #| fig-subcap: #| - "DUODENUM" #| - "JEJUNUM" #| - "ILEUM" #| - "CECUM" #| eval: true knitr::include_graphics(paste0(caminho,'Heatmap_DA_DESeq2_DUODENUM','.svg')) knitr::include_graphics(paste0(caminho,'Heatmap_DA_DESeq2_JEJUNUM','.svg')) knitr::include_graphics(paste0(caminho,'Heatmap_DA_DESeq2_ILEUM','.svg')) knitr::include_graphics(paste0(caminho,'Heatmap_DA_DESeq2_CECUM','.svg')) ``` ```{r, label="aheatmap"} #| echo: false #| warnings: false #| include: false #| eval: false #see http://www.cookbook-r.com/Graphs/Colors_%28ggplot2%29/ cols <- brewer.pal(4, "BuGn") # sets how many colours of the palette are fixed?? # can not be because it allows later even pla(2) so a number smaller than four pal <- colorRampPalette(cols) # pal is now a function colourmap <- pal(20) # makes 20 colours from the function pal counts.full <- as(phyloseq::otu_table(project.rarefied.genus.ps, taxa_are_rows=TRUE), "matrix") lineage <- apply((as(phyloseq::tax_table(project.rarefied.genus.ps),"matrix")), 1, function(x) { return(paste(x[c('Kingdom','Phylum','Class','Order','Family','Genus')],collapse=";")) } ) rownames(counts.full) <- as.character(lineage[rownames(counts.full)]) smallestGroupSize <- 1 smallestLibSize <- 10 smallestCount <- 10 keep_rows <- rowSums(counts.full >= smallestCount) >= smallestGroupSize keep_cols <- colSums(counts.full) >= smallestLibSize counts.vsd <- varianceStabilizingTransformation(counts.full[keep_rows,keep_cols],blind=TRUE) for (s in segment.item) { cols <- rownames(metadata[which(metadata$SEGMENT==s),]) pass.cols <- intersect(cols,colnames(counts.vsd)) nonzero <- apply(counts.vsd[,pass.cols],1,function(x) { sum(x)!=0 }) sample.annotation = data.frame(SEGMENT=as.factor(as.data.frame(sample_data(project.ps)[pass.cols,'SEGMENT'])$SEGMENT), TREATMENT=as.factor(as.data.frame(sample_data(project.ps)[pass.cols,'TREATMENT'])$TREATMENT), SEX=as.factor(as.data.frame(sample_data(project.ps)[pass.cols,'SEX'])$SEX) ) svg(filename=paste0(caminho,'aHeatmap_',s,'.svg'),width = 20,height=30 ) aheatmap(data.matrix(counts.vsd[nonzero,pass.cols]), col=colourmap, distfun=function(c) {vegdist(c, method="euclidean")}, hclustfun=function(x) {hclust(x, method="ward.D2")}, annCol=sample.annotation, annColors = list(SEGMENT=setNames(segment.item.color, as.character(segment.item))[levels(sample.annotation$SEGMENT)], TREATMENT=setNames(treatment.item.color, as.character(treatment.item))[levels(sample.annotation$TREATMENT)], SEX=setNames(sex.item.color, as.character(sex.item))[levels(sample.annotation$SEX)] ) ) graphics.off() } ``` Outros heatmaps foram gerados com anotações das amostras (grupos) e com todas os gêneros que tiveram abundância normalizada e possuem mais do que zero de contagem nas amostras de cada tempo (@fig-aheatmap). Os dados de abundância das amostras exibidos nesses heatmaps foram previamente transformados utlizando a transformação VST (*variance stabilizing transformation*) desconsiderando o desenho experimental (função *varianceStabilizingTransformation()* do pacote DESeq2 com parâmetro *blind=TRUE*). Os dendrogramas foram construídos a partir do algoritmo de agrupamento hierárquico com distância *euclideana* e método de ligação *ward.D2*. ```{r, label="aheatmap-shown"} #| echo: false #| eval: true #| warning: false #| lightbox: true #| layout-ncol: 1 #| layout-nrow: 4 #| label: fig-aheatmap #| fig-cap: Heatmap of genus abundances among samples by SEGMENT #| fig-subcap: #| - "DUODENUM" #| - "JEJUNUM" #| - "ILEUM" #| - "CECUM" knitr::include_graphics(paste0(caminho,'aHeatmap_DUODENUM','.svg')) knitr::include_graphics(paste0(caminho,'aHeatmap_JEJUNUM','.svg')) knitr::include_graphics(paste0(caminho,'aHeatmap_ILEUM','.svg')) knitr::include_graphics(paste0(caminho,'aHeatmap_CECUM','.svg')) ``` ```{r, label="save-protocol"} #| echo: false #| warnings: false #| eval: false save.image(file = paste0(caminho,'Analysis_Amplicon_RDP.RData')) ```