Children faecal microbiota

Author

Daniel Guariz Pinheiro

Amostras e objetivos

Este estudo tem como objetivo compreender o microbioma da microbiota intestinal de crianças por meio de amostras de fezes.

Nomenclatura

As amostras do segundo sequenciamento foram nomeadas de acordo com o padrão:

CHILD<SEQUENCING><BIOLOGICAL REPLICATE>

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

Resultados da Avaliação de qualidade (prévia)

MultiQC Report - R1 MultiQC Report - R2

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.

Os primers utilizados na verificação foram:

DNAStringSet object of length 2:
    width seq
[1]    17 CCTACGGGNGGCWGCAG
[2]    21 GACTACHVGGGTATCTAATCC

Resultado da verificação dos iniciadores (primers)

Arquivos

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.

Resultados das estimativas de tamanhos e qualidades das leituras

Arquivos *_trimSizes.txt

Arquivos *_eestats.txt

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.

Resultado das estimativas das médias de erros esperados (reads)

Arquivo 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.

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.

Resultados da avaliação de qualidade das reads após a poda de qualidade

MultiQC Report - R1 MultiQC Report - R2

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.

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.

DNAStringSet object of length 2:
    width seq
[1]    17 CCTACGGGNGGCWGCAG
[2]    17 ATTACCGCGGCTGCTGG

Quando há mais do que um amplicon por leitura, então a maior é a selecionada com um in-house Perl script GetLargestAmp.pl.

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.

Resultados da avaliação de qualidade dos amplicons

MultiQC Report - Amplicons

Estimativas de tamanhos e qualidades das leituras (amplicons)

Procedimento realizado exatamente como descrito anteriormente para as reads.

Resultados das estimativas de tamanhos e qualidades das leituras

Arquivos *_trimSizes.txt

Arquivos *_eestats.txt

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.

Resultado das estimativas das médias de erros esperados (amplicons)

Arquivo meanEE.txt

Contagens preliminares

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.

Avaliação do perfil de qualidade das leituras

As leituras R1 foram carregadas no pacote DADA2 para o reconhecimento dos amplicons. Na Figure 1 a seguir temos uma análise da qualidade das leituras em função do tamanho.

Figure 1: Quality Profile previous DADA2 processing (only one SAMPLE)

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.

O resultado (Table 1) pode também ser visualizado na Figure 2 a seguir.

Table 1: Read after filtering and trimming process

	reads.in	reads.out
CHILD21	35680	6312
CHILD22	19426	1773
CHILD24	85777	26562
CHILD25	99140	34420

Figure 2: Quality Profile post DADA2 filtering and trimming (only one sample)

Quality Profile post DADA2 filtering and trimming (full image)

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 Figure 3. 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).

Figure 3: 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.

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.

Identificação e filtragem de ASVs quiméricas

O resultado foi de 397 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.

Foram identificadas 2 ASVs quiméricas.

Números de leituras iniciais e abundâncias de ASVs

Abaixo (Table 2) os números de leituras iniciais e de abundâncias de ASVs após os processamentos até o momento.

Table 2: Reads and ASV abundances

ID	Raw.reads	Cleaned.reads	FilteredAndTrimmed.reads	Denoised.asvs.abundance	ChimeraFree.asvs.abundance	ContaminantFree.asvs.abundance	PercentUsable.reads
CHILD21	38631	38594	6312	6207	6207	6207	16.07
CHILD22	20530	20516	1773	1690	1690	1690	8.23
CHILD24	91535	91506	26562	26180	26164	26164	28.58
CHILD25	103848	103809	34420	34077	34072	34072	32.81

Tabela de contagem ASVs

O número total de ASVs válidas é de 395. As sequências fasta das ASVs podem ser acessadas no link abaixo.

Download 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%.

Identificação e filtragem de ASVs contaminantes

Além disso foram feitas filtragens de ASVs contaminantes ou não classificadas (Padrões buscados: “Chloroplast;”,“Mitochondria;”,“^Unclassified;”,“^Eukaryota;”) resultando em 395 ASVs válidas livres de contaminantes.

As sequências fasta livre de contaminantes podem ser acessadas no link abaixo.

Download ContaminantFree_ASVs.fa

Resultados das atribuições taxonômicas

As Table 3 e Table 4 sumarizam o resultado das atribuições taxonômicas correspondentes aos diferentes bancos de dados considerados no processo.

Table 3: Somatório da abundância de ASVs com atribuição até determinado nível taxonômico

	Database	Kingdom	Phylum	Class	Order	Family	Genus	Species
2	GTDB	395	374	372	372	363	335	16
3	RDP	395	385	368	368	332	254	395
4	RefRDP	395	383	364	364	327	314	9
5	Silva	395	391	391	389	358	335	395

Table 4: ASVs com atribuição até determinado nível taxonômico

	Database	UnqKingdom	UnqPhylum	UnqClass	UnqOrder	UnqFamily	UnqGenus	UnqSpecies
2	GTDB	1	7	7	17	28	48	16
3	RDP	1	4	9	13	25	43	18
4	RefRDP	1	4	9	13	24	39	9
5	Silva	1	5	7	15	26	39	17

Planilha das ASVs com atribuição taxonômica

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.

Download master.xlsx

No geral, houve 82.74% de dados com identificação da família, 47.65% de dados com informação de gênero, e 3.5% de dados com informação de espécie.

A tabela anterior (Table 2) com os números de leituras portanto foi atualizada (Table 5) com novas colunas.

Table 5: Reads and ASVs abundances (updated)

ID	Raw.reads	Cleaned.reads	FilteredAndTrimmed.reads	Denoised.asvs.abundance	ChimeraFree.asvs.abundance	ContaminantFree.asvs.abundance	PercentUsable.reads
CHILD21	38631	38594	6312	6207	6207	6207	16.07
CHILD22	20530	20516	1773	1690	1690	1690	8.23
CHILD24	91535	91506	26562	26180	26164	26164	28.58
CHILD25	103848	103809	34420	34077	34072	34072	32.81

Análise Filogenética

As 395 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 (Callahan et al. 2016) exceto pelo método de rearranjo dos ramos.

Arquivos para exploração com MicrobiomeAnalyst

Os arquivos a seguir servem para análises em plataformas externas como o MicrobiomeAnalyst. 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.

Download MA_counts.txt

Download MA_metadata.txt

Download MA_taxonomy.txt

Download MA_tree.nwk

Análises ecológicas

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 ¹ usando a função summarize_phyloseq() que sumariza os resultados do objeto phyloseq sem a rarefação (Sparsity = 0.716455696202532).

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 (5077) (Figure 4). 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 (6) foram removidas. Também podem ser consultados mais detalhes na análise de rarefação por grupos de amostras (Figure 5). As figuras foram geradas a partir da suavização dos valores de contagem de reads e respectivos intervalos de confiança (95%) quando apresentados.

(a) Rarefaction curve (original counts) by GROUP

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 (Sparsity = 0.715938303341902).

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.

Download Shared-ASVs.xlsx

(a) Percentage of shared ASVs (original dataset) by GROUP

Aglomeração de ASVs da mesma taxonomia

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 56 registros.

Análise de intersecção

As análises de intersecção a seguir estão comparando os dados das crianças (Figure 7).

Composição taxonômica

A seguir temos a composição dos mais abundantes filos e gêneros em cada amostras (Figure 8). 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.

Alfa diversidade

Os índices de alfa diversidade (“observed”, “chao1”, “diversity_inverse_simpson”, “diversity_gini_simpson”, “diversity_shannon”, “diversity_fisher”, “diversity_coverage”, “evenness_camargo”, “evenness_pielou”, “evenness_simpson”, “evenness_evar”, “evenness_bulla”, “dominance_dbp”, “dominance_dmn”, “dominance_absolute”, “dominance_relative”, “dominance_simpson”, “dominance_core_abundance”, “dominance_gini”, “rarity_log_modulo_skewness”, “rarity_low_abundance”, “rarity_rare_abundance”, “PD”) 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 (Table 6) foi incorporado também à planilha “Alpha Diversity Indices” do Excel (Section 2.13.8.3).

Para informações sobre os índices de diversidade, incluindo a forma como interpretar os resultados, consulte a seção “Indices of diversity and eveness” da página “Analysis of community ecology data in R” 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;

Download diversity.xlsx

Table 6: Alpha diversity indices

	observed	chao1	diversity_inverse_simpson	diversity_gini_simpson	diversity_shannon	diversity_fisher	diversity_coverage	evenness_camargo	evenness_pielou	evenness_simpson	evenness_evar	evenness_bulla	dominance_dbp	dominance_dmn	dominance_absolute	dominance_relative	dominance_simpson	dominance_core_abundance	dominance_gini	rarity_log_modulo_skewness	rarity_low_abundance	rarity_rare_abundance	PD
CHILD21	69	69.0	41.53	0.9759	3.918	11.292	16	0.7816	0.9255	0.6018	0.3991	0.6922	0.0620	0.1062	315	0.0620	0.0241	0.0955	0.8988	0.7442	0.0081	0.0081	62.91
CHILD24	175	175.8	95.37	0.9895	4.831	35.140	37	0.4892	0.9354	0.5450	0.5097	0.6820	0.0335	0.0648	170	0.0335	0.0105	0.0244	0.7481	1.2571	0.0410	0.0366	71.42
CHILD25	161	161.5	85.58	0.9883	4.730	31.673	32	0.4742	0.9309	0.5316	0.5408	0.6606	0.0368	0.0670	187	0.0368	0.0117	0.0416	0.7746	1.3042	0.0315	0.0309	58.71
CHILD22	37	37.0	20.27	0.9507	3.225	6.682	8	0.9119	0.8933	0.5478	0.4509	0.6218	0.0911	0.1817	154	0.0911	0.0493	0.1396	0.9509	1.1575	0.0018	0.0018	48.93

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 Figure 9 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 (Figure 10), ou seja, são os mesmos índices de diversidade alfa porém considerando todo o conjunto de amostras consideradas no sequenciamento.

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 (Figure 11). Para os perfis de diversidade utilizamos uma árvore ultramétrica a partir da coerção da árvore filogenética prévia.

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 e métodos de ordenação neste link. A diversidade beta observada pode ser avaliada por meio das figuras a seguir (Figure 12).

(a) Bray-curtis distance colour by GROUP

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).

Essa comparação foi realizada utilizando como entrada a aglomeração dos ASVs em gêneros e foram filtrados 17 gêneros, sem que os grupos de amostras tivessem a soma das abundâncias mínima de 10. Assim como bibliotecas sem um mínimo de soma de abundâncias (10), também foram removidas (). 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 (Figure 13). 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:

Download Differential_Abundance.xlsx

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 (Figure 14). 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.

References

Callahan, Ben J., Kris Sankaran, Julia A. Fukuyama, Paul J. McMurdie, and Susan P. Holmes. 2016. “Bioconductor Workflow for Microbiome Data Analysis: From Raw Reads to Community Analyses.” F1000Research 5 (November): 1492. https://doi.org/10.12688/f1000research.8986.2.

Footnotes

esparsa: valores espalhados, ou seja, com muitos zeros na matriz de abundâncias.↩︎