--- title: "Children faecal 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) obj.tmp[which(is.na(obj.tmp$down)),'down'] <- obj.tmp[which(is.na(obj.tmp$down)),'value'] obj.tmp[which(is.na(obj.tmp$up)),'up'] <- obj.tmp[which(is.na(obj.tmp$up)),'value'] # 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 = as.integer(l.sample_i), j = as.integer(c.sample_i), value = (n.l_c/n.l)*100) data.table::set( x = dt.perc.obj, i = as.integer(c.sample_i), j = as.integer(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 da microbiota intestinal de crianças por meio de amostras de fezes. ```{r, label="setup-R", include=FALSE} #| echo: false #| warnings: false #| eval: true # Primeiros passos base.dir <- "/data/MicrobiomaAraraquara/Children" Sys.setenv(BASE_DIR=base.dir) primers.file <- paste0(base.dir,"/refs/primers.fa") primers <- readFasta(primers.file) Sys.setenv(PRIMERS_FILE=primers.file) ``` ### Nomenclatura As amostras do segundo sequenciamento foram nomeadas de acordo com o padrão: CHILD\\ ## 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 \ -o ${outdir}/${bn}_cleaned_R1.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_R1.out.log \ 2>${outdir}/${bn}_cleaned_R1.err.log fastp -i ${indir}/${bn}_R2.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}_R2.html \ --json ${outdir}/${bn}_R2.json \ > ${outdir}/${bn}_cleaned_R2.out.log \ 2>${outdir}/${bn}_cleaned_R2.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) ::: {#warn-R2 .callout-warning} #### Decisão As leituras R2 do sequenciamento 2 foram praticamente eliminadas e a partir deste momento deixarão de fazer parte das análises. Os *amplicons* a partir daqui serão representados apenas pelas *reads* R1. ::: ```{bash, label='flash-merge-fake', 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}/*_R1*.fastq) ; do bn=$(basename ${fq} _cleaned_R1.fastq) echo "Linking processed R1 file with amplicon file ..." abssrc=`readlink -f ${fq}` absdest=`readlink -f ${outdir}/${bn}.fastq` ln -s ${abssrc} ${absdest} 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. Porém, como somente foram consideradas as *reads* R1, também consideramos apenas o *primer* *forward* e o *primer* *reverse* foi considerado além do oligo CCTACGGGNGGCWGCAG (341F), também o oligo ATTACCGCGGCTGCTGG (534R), que permitem a captura da região hipervariável V3. Um máximo de 5 diferenças foi considerado como uma correspondência (*match*) aceito entre as sequências dos *primers* e as sequências das *reads*. ```{r, label="print-newprimers"} #| echo: false #| eval: true print(newprimers <- c(primers[1]@sread,DNAStringSet("ATTACCGCGGCTGCTGG"))) ``` ```{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} echo ${PRIMERS_FILE} fwdprimer=`grep '^>' -m 1 -A 1 ${PRIMERS_FILE} | tail -1` #revprimer=`grep '^>' -m 2 -A 1 ${PRIMERS_FILE} | tail -1` #V3 revprimer="ATTACCGCGGCTGCTGG" new_primers_file="${outdir}/primerdb.fa" echo ">FWD" > ${new_primers_file} echo "${fwdprimer}" >> ${new_primers_file} echo ">REV" >> ${new_primers_file} echo "${revprimer}" >> ${new_primers_file} minamp=100 maxamp=300 for fq in `ls ${indir}/*.fastq 2>/dev/null`; do bn=`basename ${fq} .fastq`; usearch11 -search_pcr ${fq} \ -db ${new_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 um *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\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}') vld=$(cat ${outdir}/amplicons/${bn}.fastq | wc -l | sed 's/ /\t/g' | awk '{print $1/4}') echo -e "${id}\t${raw}\t${cln}\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} #| echo: false #| eval: false # Criando um atalho para o caminho que iremos salvar as saídas: seqn<-NULL caminho <- paste0(base.dir, "/output/dada2/",ifelse(!is.null(seqn),seqn,'')) 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 input.sel <- input[grep(paste0('^CHILD',ifelse(!is.null(seqn),seqn,'')),input$nomes),] # Carregando os caminhos dos arquivos reads <- input.sel$reads # Separando nome das amostras: nomes <- input.sel$nomes # Criando caminhos de saída para as bibliotecas processadas com dada2: names(reads) <- nomes 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 if (file.exists(paste0(caminho,'/Analysis_Amplicon_RDP_4.RData'))) { load(file = paste0(caminho,'/Analysis_Amplicon_RDP_4.RData')) } if (is.null(rand.samp)) { rand.samp <- nomes[ sample(1:length(reads),1) ] } ``` #### Avaliação do perfil de qualidade das leituras 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')) ``` #### Filtragem e poda de leituras A seguir foi realizada uma filtragem e poda nas sequências de acordo com o perfil de qualidade observado, dessa forma, consideramos truncar as sequências até um determinado tamanho (*truncLen = 180*), com um valor máximo de expectativa de erros (*maxEE*) de 1, 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 = 1, maxN = 0, rm.phix = T, compress = F, multithread=T, rm.lowcomplex=5, verbose=T ) rownames(qc) <- gsub('.fastq','',rownames(qc)) reads_filt <- reads_filt[rownames(qc[which(qc[,'reads.out']>0),])] nomes <- rownames(qc[which(qc[,'reads.out']>0),]) ``` O resultado (@tbl-filttrim) pode também ser visualizado na @fig-plotQPrf a seguir. ```{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). ```{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[which(qc[,'reads.out']>0),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... if (file.exists(paste0(caminho,'master.xlsx'))) { 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, 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 #| eval: false # Salvando tabelas para o MicrobiomeAnalyst: tb_meta <- read.delim(file = paste0(base.dir,"/metadata_children.tsv"), header = T, sep = "\t", check.names=FALSE) %>% filter(`#NAME` %in% nomes) 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_data <- 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") for (n in nomes) { children.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_data)), phyloseq::sample_data(tb_meta.ma[c(paste0(n,'.1'),paste0(n,'.2'),paste0(n,'.3')),]), phyloseq::phy_tree(fitGTR$tree) ) # rarefy with replacement children.rarefied.ma = suppressMessages( rarefy_even_depth( children.ma, rngseed=1, replace=T) ) for (asv in rownames( as.data.frame(otu_table(children.rarefied.ma)) ) ) { tb_cont.ma[asv, paste0(n,'.1')] <- as.data.frame(otu_table(children.rarefied.ma))[asv, paste0(n,'.1')] tb_cont.ma[asv, paste0(n,'.2')] <- as.data.frame(otu_table(children.rarefied.ma))[asv, paste0(n,'.2')] tb_cont.ma[asv, paste0(n,'.3')] <- as.data.frame(otu_table(children.rarefied.ma))[asv, paste0(n,'.3')] } } #colSums(tb_cont.ma[,c('CHILD21.1','CHILD21.2','CHILD21.3')]) #sum(tb_cont[, 'CHILD21']) #head(tb_cont.ma[,c('CHILD21.1','CHILD21.2','CHILD21.3')]) 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 set.seed(1) metadata$GROUP <- gsub('\\d+','',metadata$SOCIOECOSTAT) project.ps <- phyloseq::phyloseq(phyloseq::otu_table(as.matrix(otu),taxa_are_rows=TRUE), phyloseq::tax_table(as.matrix(taxonomy_data)), phyloseq::sample_data(metadata), phyloseq::phy_tree(fitGTR$tree) ) 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_data)), phyloseq::sample_data( metadata[names( sample_sums( project.ps )[sample_sums(project.ps) < rarefy.threshold] ), ] ), phy_tree(project.ps) ) # rarefy without replacement project.rarefied.ps <- rarefy_even_depth( project.ps, sample.size=rarefy.threshold, replace=F, trimOTUs=F, rngseed=1, verbose=TRUE) project.final.ps <- merge_phyloseq(project.rarefied.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.rarefied.ps) summ.final.ps<-summarize_phyloseq(project.final.ps) #sparsity #(sum(otu_table(project.rarefied.ps)==0)/(dim(otu_table(project.rarefied.ps))[1]*dim(otu_table(project.rarefied.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 group.item <- c("A","B","C") group.item.color <- c("darkred","#BA8E23","darkgreen") group.item.linetype <- c("solid","solid","solid") child.item <- c("1","2","4","5") child.item.color <- c("blue","red","orange","purple") child.item.linetype <- c("solid","solid","solid","solid") ######## ### FULL ######## rareres <- ggrarecurve.get_rarecurve(project.ps, step=100, file=paste0(caminho,'Rarecurve_DADA2','.svg') ) svg(filename=paste0(caminho,'Rarecurve_GROUP_DADA2','.svg') ) prare.group <- ggrarecurve.rarecurve(obj=rareres, factorNames="GROUP", indexNames=c("Observe"), shadow=TRUE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(group.item.color, group.item))+ scale_fill_manual(values=setNames(group.item.color, group.item))+ scale_linetype_manual(values=setNames(group.item.linetype, group.item))+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="white"), strip.text.x = element_text(face="bold")) prare.group graphics.off() svg(filename=paste0(caminho,'Rarecurve_GROUP_SAMPLE_DADA2','.svg') ) prare.group.sample <- ggrarecurve.rarecurve(obj=rareres, factorNames="GROUP", indexNames=c("Observe"), shadow=FALSE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(group.item.color, group.item))+ scale_fill_manual(values=setNames(group.item.color, group.item))+ scale_linetype_manual(values=setNames(group.item.linetype, group.item))+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="white"), strip.text.x = element_text(face="bold")) prare.group.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_GROUP_DADA2_rarefied','.svg') ) prare.group.rarefied <- ggrarecurve.rarecurve(obj=rareres.rarefied, factorNames="GROUP", indexNames=c("Observe"), shadow=TRUE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(group.item.color, group.item))+ scale_fill_manual(values=setNames(group.item.color, group.item))+ scale_linetype_manual(values=setNames(group.item.linetype, group.item))+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="white"), strip.text.x = element_text(face="bold")) prare.group.rarefied graphics.off() svg(filename=paste0(caminho,'Rarecurve_GROUP_SAMPLE_DADA2_rarefied','.svg') ) prare.group.sample.rarefied <- ggrarecurve.rarecurve(obj=rareres.rarefied, factorNames="GROUP", indexNames=c("Observe"), shadow=FALSE, se=FALSE, linesize=1 ) + scale_color_manual(values=setNames(group.item.color, group.item))+ scale_fill_manual(values=setNames(group.item.color, group.item))+ scale_linetype_manual(values=setNames(group.item.linetype, group.item))+ theme_bw()+ theme(axis.text=element_text(size=8), panel.grid=element_blank(), strip.background = element_rect(colour=NA,fill="white"), strip.text.x = element_text(face="bold")) prare.group.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 #| eval: true #| warning: false #| lightbox: true #| layout-ncol: 2 #| layout-nrow: 2 #| label: fig-rarecurve-groups #| fig-cap: Rarefaction analysis by groups #| fig-subcap: #| - "Rarefaction curve (original counts) by GROUP" #| - "Rarefaction curve (rarefied counts) by GROUP" #| - "Rarefaction curve (original counts) by SAMPLE" #| - "Rarefaction curve (rarefied counts) by SAMPLE" knitr::include_graphics(paste0(caminho,'Rarecurve_GROUP_DADA2','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_GROUP_DADA2_rarefied','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_GROUP_SAMPLE_DADA2','.svg')) knitr::include_graphics(paste0(caminho,'Rarecurve_GROUP_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(GROUP=as.factor(sample_data(project.ps)[['GROUP']])) 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(GROUP=setNames(group.item.color, as.character(group.item))[levels(otu.perc.sample.annotation$GROUP)]) ) 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(GROUP=setNames(group.item.color, as.character(group.item))[levels(otu.perc.sample.annotation$GROUP)]) ) graphics.off() ``` ```{r, label="shared-otu-spreadsheet"} #| echo: false #| warnings: false #| eval: true xfun::embed_file(paste0(caminho, "Shared-ASVs.xlsx")) ``` ```{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 = phyloseq::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 após procedimento de rarefação temos `r dim(otu_table(project.rarefied.genus.ps))[1]` registros. #### Análise de intersecção #### Análise de intersecção As análises de intersecção a seguir estão comparando os dados das crianças (@fig-upset-by-child). ```{r, label="upset-plots"} #| echo: false #| warning: false #| include: false #| eval: false metadata_factors<-sample_data(project.ps) metadata_factors$GROUP <- as.factor(as.character(metadata_factors$GROUP)) metadata_factors$CHILD <- as.factor(as.character(metadata_factors$CHILD)) metadata.grp <- metadata_factors upsetda <- get_upset(obj=as.data.frame(t(as.matrix(otu_table(project.final.ps)[,rownames(metadata.grp)]))), sampleda=metadata.grp, factorNames='CHILD') 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"))[1:length(unique(comb_degree(m1)))], 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.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.grp)]))), sampleda=metadata.grp, factorNames='CHILD') 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"))[1:length(unique(comb_degree(m2)))], 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.svg') ) plot(upset.genus) graphics.off() ``` ```{r, label="upset-plots"} #| echo: false #| eval: true #| lightbox: true #| layout-ncol: 2 #| layout-nrow: 1 #| label: fig-upset-by-child #| fig-cap: Intersection analyses (UpSet plots) among GROUPs #| fig-subcap: #| - "ASV level" #| - "Genus level" knitr::include_graphics(paste0(caminho,'UpSet_asvs','.svg')) knitr::include_graphics(paste0(caminho,'UpSet_genus','.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) #plot_bar(project.final.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("GROUP"), 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=25 ) phylumPlot graphics.off() compGenus <- trans_abund$new(dataset = meco, taxrank = "Genus", ntaxa = 10) genusPlot <- compGenus$plot_bar(others_color = "grey90", facet = c("GROUP"), 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=25 ) genusPlot graphics.off() ``` ```{r, label="taxonomy-composition-shown"} #| echo: false #| eval: true #| warning: false #| lightbox: true #| layout-ncol: 2 #| label: fig-compo #| fig-cap: Compositional plot #| fig-subcap: #| - "Phylum level (Top 5)" #| - "Genus level (Top 10)" 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 #| include: 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") 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 = 50/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. Considerando cada amostra um ambiente, podemos considerar que a diversidade relacionada ao conjunto das amostras (CRIANÇAS) 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"} #| eval: false #| echo: false #| warning: false #| include: false metadata.div <- metadata_factors for (i in colnames(div.table)) { metadata.div[[i]] <- div.table[[i]] } for (i in colnames(div.table)) { p1 <- ggviolin(metadata.div, x = "GROUP", y = i, add = "boxplot", fill = "GROUP") + scale_fill_manual(values=setNames(group.item.color, group.item), breaks=group.item) p1 <- p1 + geom_pwc(method="wilcox_test", label = "{p.format} {p.signif}", hide.ns = 'p') + scale_x_discrete(guide = guide_axis(angle = 45)) ggsave(paste0(caminho,'AlphaDivPlot_GROUP_CHILD_',i,'.svg'), p1, width=20, unit='cm' ) p2 <- ggviolin(metadata.div, x = "CHILD", y = i, add = "boxplot", fill = "GROUP") p2 <- p2 + geom_pwc(method="wilcox_test", label = "{p.format} {p.signif}", hide.ns = 'p') + scale_x_discrete(guide = guide_axis(angle = 45)) ggsave(paste0(caminho,'AlphaDivPlot_CHILD_GROUP_',i,'.svg'), p2, width=20, unit='cm' ) p3 <- ggviolin(metadata.div, x = "SEQUENCING", y = i, add = "boxplot", fill = "gray") ggsave(paste0(caminho,'GammaDivPlot_',i,'.svg'), p3, width=20, unit='cm' ) } ``` ```{r} #| eval: true #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 2 #| label: fig-plotalphadiv #| fig-cap: Alpha diversity indices per child #| fig-subcap: #| - "Richness index (observed)" #| - "Diversity index (shannon)" #for (i in colnames(div.table)) { #knitr::include_graphics(paste0(caminho,'AlphaDivPlot_',i,'.svg')) #} knitr::include_graphics(paste0(caminho,'AlphaDivPlot_GROUP_CHILD_observed','.svg')) knitr::include_graphics(paste0(caminho,'AlphaDivPlot_GROUP_CHILD_diversity_shannon','.svg')) ``` ```{r} #| eval: true #| 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)" #for (i in colnames(div.table)) { #knitr::include_graphics(paste0(caminho,'GammaDivPlot_',i,'.svg')) #} 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). Para os perfis de diversidade utilizamos uma árvore ultramétrica a partir da coerção da árvore filogenética prévia. ```{r, label="hill-profile"} #| echo: false #| warning: false #| include: false #| eval: false ultrametric_tree <- force.ultrametric( phy_tree(project.final.ps) ) group <- metadata[rownames(sample_data(project.final.ps)),'GROUP'] samp <- rownames(sample_data(project.final.ps)) hierarchy.group <- data.frame(Sample=samp, Group=group) svg(paste0(caminho,'Hill_profile_alpha_GROUP','.svg') ) tmp_table <- as.data.frame(otu_table(project.final.ps)) few.reps <- hierarchy.group[which( hierarchy.group$Group %in% names(table(hierarchy.group$Group)[table(hierarchy.group$Group)<2]) ),'Sample'] tmp_table <- as.data.frame(otu_table(project.final.ps)[,setdiff(colnames(tmp_table),few.reps)]) hierarchy.group <- subset(hierarchy.group, ! Sample %in% few.reps ) for (fr in few.reps) { tmp_table[rownames(otu_table(project.final.ps)) ,paste0(fr,'.1')] <- as.numeric( otu_table(project.final.ps)[,fr] ) tmp_table[rownames(otu_table(project.final.ps)) ,paste0(fr,'.2')] <- as.numeric( otu_table(project.final.ps)[,fr] ) hierarchy.group <- rbind(hierarchy.group,c(Sample=paste0(fr,'.1'),Group=metadata[fr,'GROUP'])) hierarchy.group <- rbind(hierarchy.group,c(Sample=paste0(fr,'.2'),Group=metadata[fr,'GROUP'])) } div_profile_plot(div_profile(count=tmp_table, hierarchy=hierarchy.group ),colour=setNames(group.item.color, group.item)) rm(tmp_table) graphics.off() # BY CHILD hierarchy.child <- data.frame(Sample=rownames(sample_data(project.final.ps)), Group=paste0("CHILD",metadata[rownames(sample_data(project.final.ps)),'CHILD'])) svg(paste0(caminho,'Hill_profile_alpha_CHILD','.svg') ) tmp_table <- as.data.frame(otu_table(project.final.ps)) few.reps <- hierarchy.child[which( hierarchy.child$Group %in% names(table(hierarchy.child$Group)[table(hierarchy.child$Group)<2]) ),'Sample'] if (length(colnames(tmp_table))==length(few.reps)) { tmp_table <- data.frame() } else { tmp_table <- as.data.frame(otu_table(project.final.ps)[,setdiff(colnames(tmp_table),few.reps)]) } hierarchy.child <- subset(hierarchy.child, ! Sample %in% few.reps ) for (fr in few.reps) { tmp_table[rownames(otu_table(project.final.ps)) ,paste0(fr,'.1')] <- as.numeric( otu_table(project.final.ps)[,fr] ) tmp_table[rownames(otu_table(project.final.ps)) ,paste0(fr,'.2')] <- as.numeric( otu_table(project.final.ps)[,fr] ) hierarchy.child <- rbind(hierarchy.child,c(Sample=paste0(fr,'.1'),Group=paste0("CHILD",metadata[fr,'CHILD']))) hierarchy.child <- rbind(hierarchy.child,c(Sample=paste0(fr,'.2'),Group=paste0("CHILD",metadata[fr,'CHILD']))) } colnames(hierarchy.child) <- c('Sample','Group') div_profile_plot(div_profile(count=tmp_table, hierarchy=hierarchy.child ),colour=setNames(child.item.color, paste0("CHILD",child.item))) rm(tmp_table) 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"} #| eval: true #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 3 #| label: fig-hillprofile #| fig-cap: Hill diversity plots #| fig-subcap: #| - "Alpha diversity profiles by GROUP" #| - "Alpha diversity profiles by CHILD" #| - "Gamma diversity profile" knitr::include_graphics(paste0(caminho,'Hill_profile_alpha_GROUP','.svg')) knitr::include_graphics(paste0(caminho,'Hill_profile_alpha_CHILD','.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). ```{r, label="beta-div-plot"} #| echo: false #| warning: false #| include: false #| eval: false meco$cal_betadiv(unifrac = TRUE) t.beta.group.bray <- trans_beta$new(dataset = meco, group = "GROUP", measure = "bray") t.beta.group.bray$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_bray_GROUP','.svg'), t.beta.group.bray$plot_ordination(plot_color="GROUP", color_values = setNames(group.item.color, group.item), plot_type = c("point")), width=20, unit='cm' ) t.beta.group.jaccard <- trans_beta$new(dataset = meco, group = "GROUP", measure = "jaccard") t.beta.group.jaccard$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_jaccard_GROUP','.svg'), t.beta.group.jaccard$plot_ordination(plot_color="GROUP", color_values = setNames(group.item.color, group.item), plot_type = c("point")), width=20, unit='cm' ) t.beta.group.wei_unifrac <- trans_beta$new(dataset = meco, group = "GROUP", measure = "wei_unifrac") t.beta.group.wei_unifrac$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_wei_unifrac_GROUP','.svg'), t.beta.group.wei_unifrac$plot_ordination(plot_color="GROUP", color_values = setNames(group.item.color, group.item), plot_type = c("point")), width=20, unit='cm' ) t.beta.group.unwei_unifrac <- trans_beta$new(dataset = meco, group = "GROUP", measure = "unwei_unifrac") t.beta.group.unwei_unifrac$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_unwei_unifrac_GROUP','.svg'), t.beta.group.unwei_unifrac$plot_ordination(plot_color="GROUP", color_values = setNames(group.item.color, group.item), plot_type = c("point")), width=20, unit='cm' ) meco$sample_table$CHILD<-as.factor(meco$sample_table$CHILD) t.beta.child.bray <- trans_beta$new(dataset = meco, group = "CHILD", measure = "bray") t.beta.child.bray$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_bray_CHILD','.svg'), t.beta.child.bray$plot_ordination(plot_color="CHILD", color_values = setNames(child.item.color, child.item), plot_type = c("point")), width=20, unit='cm' ) t.beta.child.jaccard <- trans_beta$new(dataset = meco, group = "CHILD", measure = "jaccard") t.beta.child.jaccard$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_jaccard_CHILD','.svg'), t.beta.child.jaccard$plot_ordination(plot_color="CHILD", color_values = setNames(child.item.color, child.item), plot_type = c("point")), width=20, unit='cm' ) t.beta.child.wei_unifrac <- trans_beta$new(dataset = meco, group = "CHILD", measure = "wei_unifrac") t.beta.child.wei_unifrac$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_wei_unifrac_CHILD','.svg'), t.beta.child.wei_unifrac$plot_ordination(plot_color="CHILD", color_values = setNames(child.item.color, child.item), plot_type = c("point")), width=20, unit='cm' ) t.beta.child.unwei_unifrac <- trans_beta$new(dataset = meco, group = "CHILD", measure = "unwei_unifrac") t.beta.child.unwei_unifrac$cal_ordination(method = "PCoA") ggsave(paste0(caminho,'PCoA_unwei_unifrac_CHILD','.svg'), t.beta.child.unwei_unifrac$plot_ordination(plot_color="CHILD", color_values = setNames(child.item.color, child.item), plot_type = c("point")), width=20, unit='cm' ) ``` ```{r, label="beta-div-plot-shown"} #| eval: true #| echo: false #| warning: false #| lightbox: true #| layout-ncol: 4 #| layout-nrow: 2 #| label: fig-betadiv-sample #| fig-cap: Beta diversity analyses using PCoA ordination (SAMPLE) #| fig-subcap: #| - "Bray-curtis distance colour by GROUP" #| - "Jaccard distance colour by GROUP" #| - "Weighted UniFrac distance by GROUP" #| - "Unweighted UniFrac distance by GROUP" #| - "Bray-curtis distance colour by CHILD" #| - "Jaccard distance colour by CHILD" #| - "Weighted UniFrac distance by CHILD" #| - "Unweighted UniFrac distance by CHILD" knitr::include_graphics(paste0(caminho,'PCoA_bray_GROUP','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_jaccard_GROUP','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_wei_unifrac_GROUP','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_unwei_unifrac_GROUP','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_bray_CHILD','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_jaccard_CHILD','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_wei_unifrac_CHILD','.svg')) knitr::include_graphics(paste0(caminho,'PCoA_unwei_unifrac_CHILD','.svg')) ``` ### Análise de abundância diferencial Aqui podemos fazer uma comparação da abundância dos *taxa* (gêneros) entre as microbiotas nos grupos amostrais em um determinado tempo. 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 metadata.deg <- metadata_factors prefix.group <- list('A'='A', 'B'='B', 'C'='C' ) counts.deg <- as(phyloseq::otu_table(project.rarefied.genus.ps, taxa_are_rows=TRUE), "matrix")[,rownames(metadata.deg)] 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)]) dds <- DESeqDataSetFromMatrix(countData = counts.deg, colData = DataFrame( metadata.deg[colnames(counts.deg),] ), design = ~0+GROUP, tidy=FALSE) 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 <- varianceStabilizingTransformation(dds, blind=FALSE) # #rld <- rlog(dds, blind=FALSE) # #meanSdPlot(assay(ntd)) # #meanSdPlot(vsd) # #meanSdPlot(assay(rld)) counts.ntd <- as.data.frame(assay(ntd)) #counts.ntd <- as.data.frame(assay(vsd)) letras <- list() file.remove(paste0(caminho,'Differential_Abundance.xlsx')) wb <- createWorkbook() s <- 1 for (grp in levels(ddsClean$GROUP)) { smpst <- rownames(subset(metadata.deg[names(keep_cols[keep_cols]),], GROUP==grp)) if (length(smpst)>1) { counts.ntd[[paste0(grp,'.','sum')]] <- rowSums(counts.ntd[,smpst]) counts.ntd[[paste0(grp,'.','mean')]] <- rowSums(counts.ntd[,smpst])/length(smpst) } else { counts.ntd[[paste0(grp,'.','sum')]] <- counts.ntd[,smpst] counts.ntd[[paste0(grp,'.','mean')]] <- counts.ntd[,smpst]/length(smpst) } } # resultsNames(ddsClean) all.comb.s <- apply(t( combn(paste0('GROUP',levels( metadata.deg$GROUP) ),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('GROUP'),'',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('^.*GROUP','',ctrst[1]) t2 <- gsub('^.*GROUP','',ctrst[2]) metadata.deg.tmp <- metadata.deg[names(keep_cols[keep_cols]),] smps1 <- rownames(subset(metadata.deg.tmp, GROUP==t1)) smps2 <- rownames(subset(metadata.deg.tmp, GROUP==t2)) if (( length(smps1) >= 1 ) && ( length(smps2) >= 1 )) { da.res <- (merge(x=dalrt.df, y=counts.ntd[,c(smps1, smps2, paste(t1,'mean',sep='.'), smps2, paste(t2,'mean',sep='.') )], by='row.names'))[,c('Row.names', smps1, paste(t1,'mean',sep='.'), smps2, paste(t2,'mean',sep='.'), 'log2FoldChange', 'pvalue', 'padj' )] rownames(da.res) <- da.res$Row.names wb.name <- paste0("DA_",s,'_',paste(c(as.character(prefix.group[t1]),as.character(prefix.group[t2])),collapse='-')) addWorksheet(wb, wb.name) writeData(wb, wb.name, da.res, startRow = 1, startCol = 1) for (g in rownames(compara)) { compara[g,new.comb.s] <- ifelse(is.na(da.res[g,'padj']), FALSE, ifelse(da.res[g,'padj']<=0.1, 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$GROUP)] <- difcomparaL[levels(ddsClean$GROUP)] } } } 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 s<-1 if ( ! is.null(letras[[s]]) ) { cols <- paste0(levels(ddsClean$GROUP),'.mean') any.signif <- apply(letras[[s]], 1, function(x) { length(unique(x))>1 } ) sum.ordered <- order(counts.ntd$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, # clustering_distance_rows=function(c) {vegdist(c, method="bray")}, # clustering_distance_columns=function(c) {vegdist(c, method="bray")}, # clustering_method_rows="ward.D2", # clustering_method_columns="ward.D2", 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.svg'),width = 20,height=30 ) #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() } #res[[ 'Colitis-Colitis_personalized_treatment' ]][grep('Duncaniella',rownames(letras[[s]])), ] #letras[[s]][grep('Duncaniella',rownames(letras[[s]])),levels(ddsClean$GROUP)] ``` ```{r, label="heatmap-shown"} #| echo: false #| eval: true #| warning: false #| lightbox: true #| layout-ncol: 1 #| layout-nrow: 1 #| label: fig-heamapda #| fig-cap: Heatmap of significant genus abundances among sample GROUPS #| fig-subcap: #| - "Socieconomic status" knitr::include_graphics(paste0(caminho,'Heatmap_DA_DESeq2','.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) cols <- rownames(metadata_factors) pass.cols <- intersect(cols,colnames(counts.vsd)) nonzero <- apply(counts.vsd[,pass.cols],1,function(x) { sum(x)!=0 }) sample.annotation = data.frame(GROUP=metadata_factors[pass.cols,'GROUP']) 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(GROUP=setNames(group.item.color, as.character(group.item))[levels(sample.annotation$GROUP)]), filename=paste0(caminho,'aHeatmap.svg'), width=20, height=20, fontsize=12, cellwidth = 40, cellheight = 20 ) grid.text(label="Heatmap", x = unit(0.5, "npc"), y = unit(0.5, "npc"), just = "centre", hjust = NULL, vjust = NULL, rot = 0, check.overlap = FALSE, default.units = "npc", name = NULL, gp = gpar(), draw = TRUE, vp = NULL) ``` 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: 1 #| label: fig-aheatmap #| fig-cap: Heatmap of genus abundances among samples #| fig-subcap: #| - "Socioeconomic status" knitr::include_graphics(paste0(caminho,'aHeatmap','.svg')) ``` ```{r, label="save-protocol"} #| echo: false #| warnings: false #| eval: false save.image(file = paste0(caminho,'/Analysis_Amplicon_RDP_4.RData')) ```