# QIIME2 - gene 16S rRNA
Neste tutorial 🗒 você vai aprender 🤓 a processar e analisar dados de sequenciamento *high throughput* do gene 16S rRNA, utilizando um *pipeline open source* chamado 🔗[Qiime](https://qiime2.org/) *|chime|* (*Quantitative Insights Into Microbial Ecology*). Atualmente ele está na sua versão 2020.2, mas no nosso servidor temos instalada a última do ano passado (2019.10) ✅ que é mais estável.
## Colaboradores

* :female-scientist::female-technologist: MSc. Kelly J. Hidalgo Martinez
Microbióloga
Doutoranda em Genética e Biologia Molecular
Instituto de Biologia - UNICAMP
:iphone: Whastapp: [phone]
:mailbox_with_mail: Email: [e-mail]

* :male-scientist::male-technologist: Victor Borin Centurion
Biomédico
Doutorando em Genética e Biologia Molecular
Instituto de Biologia - UNICAMP
:iphone:Whastapp: [phone]
:mailbox_with_mail: Email: [e-mail]

* :male-scientist::male-technologist: Dr. Tiago Palladino Delforno
Biólogo
:mailbox_with_mail: Email: [e-mail]

---
## Requisitos
**Importante**:exclamation:
Não comece este tutorial 🗒 se você ainda não leu e praticou com anteriores, listados aqui:
1. Não preparou seu 💻 para trabalhar no servidor? então vai em 🔗[Preparativos para começar](https://) (Obrigatório)
2. Já tem uma ideia de que é a tela preta e como se navega nela? não? então vai em 🔗 [UNIX Shell - Linux](https://) (Obrigatório)
3. Você precisa conhecer o tipo de arquivos com os que vai trabalhar. Em 🔗[Tipos de arquivos ](https://) encontra essa informação (Obrigatório).
4. Você precisa saber como activar e desactivar ambientes virtuais em Conda para o uso das ferramentas presentes no nosso servidor. Se ainda não passou por aí vai 🔗[aqui](https://) (Obrigatório).
5. Este turorial 🗒 é uma continuação do parte I do tutorial 🗒 🔗[Control de qualidade e trimagem](https://). Então não dá para você começar aqui ⛔️. Além porque aí você vai aprender bem os parâmetros de qualidade neste tipo de dados. Vai lá primeiro então. Usando Qiime2 você consegue fazer o controle de qualidade também, mas sempre é bom conhecer diversos jeitos e ferramentas.
6. Para entender um pouco da estrutura do gene 16S e de barcoding 🔗[vai aqui](https://).

## Alguns tips interessantes antes de começar
Não esqueça de ter presentes estes tips 💁🏻‍♀️
* O tip mais importante é você ser sempre crítico 🤓 no que você está fazendo, **não** se trata de cópiar e colar comando 🙄 , tente entender o que você está fazendo, que tem por trás de cada processo, o que está acontecendo 🤔 . Este tutorial 🗒 é só um exemplo, mas tem muitos parâmetros para ser trocados segundo seus dados.
* Lembra da tecla [Tab] :keyboard: ? para autocompletar nomes de pastas 📁 e/ou arquivos e evitar erros de digitação.
* Com a seta para cima você pode voltar nos comandos que você já usou anteiormente.
* Lembre-se que a linha de comando é sensível a maiúsculas e minúsculas.
* Todo comando/programa tem um menu de ajuda **help**. `comando --help`
* Na linha de comando o simbolo `#`, significa que não vai ser executado, e é usado para deixar mensagens.
* Para descarregar 🔽 arquivos desde o servidor não esqueça de usar 🔗[Filezilla](https://)

---
## Vamos lá! :beginner:

No tutorial 🗒 de **Control de qualidade e trimagem** você aprendeu como ver gráficamente a qualidade das *reads* e a filtrar aquelas com baixa qualidade, para levar as melhores para os análises *downstream*. Bem, agora vamos continuar trabalhando com essas sequências de boa qualidade e sem primers para avaliar a diversidade taxonômica das amostras exemplo.

### Pasta de trabalho

Lembre-se que a sua 📁 de trabalho é `/data/usuário/16S_tutorial`

## Confira onde você está
pwd

Se não estiver em `/data/usuário/16S_tutorial`, então use os comandos aprendidos para entrar na pasta. Precisa lembrar comandos? :link:[Vai aqui!](https://)

### Ativar o qiime2 no conda
Se tiver dúvidas sobre como ativar ambientes no conda, 🔗[vai aqui](https://).

source /opt/Miniconda3/bin/activate qiime2-2019.10

## 4 Importar os arquivos `.fq` como *"artefato"* dentro do Qiime
Para importar os arquivos com as sequências já limpas e prontas para trabalhar dentro do Qiime2, é preciso criar um arquivo `.txt` onde vamos a colocar o nome das amostras, e os caminhos da 📁 onde estão os arquivos `.fq` forward e reverse. O arquivo vai ser chamado `ManifestFile.txt`.
Que é um artefato (`.qza`)? 🔗[Glosario qiime2](https://docs.qiime2.org/2020.2/glossary/). Por enquanto posso te adiantar que esse tipo de arquivo não pode ser "lido" por nenhum tipo de programa. Para conseguir visualizar o conteúdo destes arquivos é preciso transformar ele para um arquivo tipo *visualization* (`.qzv`).

# Abra o editor de texto nano 
nano ManifestFile.txt

Você pode copiar e colar o texto embaixo, mas leve em conta que, o arquivo `ManifestFile.txt` é um arquivo de texto separado por tabs, e só vai funcionar desse jeito, então verifique que ficou correto ✅

## Colunas separadas por tab
sample-id	forward-absolute-filepath	reverse-absolute-filepath
amostra1	$PWD/03.CleanData/amostra1_r1_paired.fq.gz	$PWD/03.CleanData/amostra1_r2_paired.fq.gz
amostra2	$PWD/03.CleanData/amostra2_r1_paired.fq.gz	$PWD/03.CleanData/amostra2_r2_paired.fq.gz
amostra3	$PWD/03.CleanData/amostra3_r1_paired.fq.gz	$PWD/03.CleanData/amostra3_r2_paired.fq.gz
amostra4	$PWD/03.CleanData/amostra4_r1_paired.fq.gz	$PWD/03.CleanData/amostra4_r2_paired.fq.gz

Deve ficar assim:

## Para sair
Ctrl + x
## Para salvar
S
Enter

Antes de rodar o primeiro comando no qiime2, vamos aprender algumas generalidades sobre a estrutura dos comandos do qiime2.
1) Sempre começam com a palavra qiime
2) Enseguida se coloca a ferramenta a ser usada, p.e. tools, dada2, feature-classification, etc). Para conhecer todas as possibilidades do qiime2 vai no 🔗[site](https://docs.qiime2.org/2019.10/) ou digite `qiime --help` na linha de comando.
3) A continuação é necessário colocar o comando a ser rodado dessa ferramenta e finalizar com `\`. A barra invertida significa que você pode escrever na linha de embaixo sem rodar o comando.
4) Posteriormente tem que colocar os parâmetros do comando que você quer usar. Alguns parâmetros são **obrigatórios**. Outros tem valores **default**, ou seja mesmo que você não coloque o parâmetro ele vai ser usado com o valor por *default* da ferramenta. Para saber quais são os parâmetros obrigatórios e os que tem valores default basta com digitar (comando genêrico) `qiime ferramenta comando --help` 🆘
5) Por último um dos parâmetros a colocar é o 📁 de saída, onde serão colocados todos arquivos de saída do comando
Tranquil@ agora mesmo pode parecer confuso, mas na medida que você vai usando os comando do qiime2, você vai conseguir "*pegar o jeito*" 😀

Agora sim o comando para importar os arquivos de sequências dentro do qiime2

## Primeiro crie um diretório para colocar o "artefato" gerado com todas as amostras
mkdir 04.ImportedReads
## Importar como artefato em qiime2
qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path ManifestFile.txt \
  --output-path 04.ImportedReads/reads_trimmed.qza \
  --input-format PairedEndFastqManifestPhred33V2
## Visualização
qiime demux summarize \
 --i-data 04.ImportedReads/reads_trimmed.qza \
 --o-visualization 04.ImportedReads/reads_trimmed.qzv

Use a plataforma qiime2 view para 👀 arquivos `.qzv` [Qiime2 View](https://view.qiime2.org)


Lembra um pacote de *scripts* que descarregou no tutorial anterior? chamado microbiome_helper. Vamos usar ele aqui para gerar uma tabela com o número de sequências de cada amostra.

/home/bioinfo/Documentos/microbiome_helper/qiime2_fastq_lengths.py 04.ImportedReads/reads_trimmed.qza --proc 4 -o read_counts.tsv
## Ver o arquivo
nano read_counts.tsv
## Complete seu arquivo QC_controle.xlsx!!

---
**Uma pausa nos comando para entender uns conceitos muito importantes!!**

Vamos revisar alguns conceitos importantes...
**OTU *(Operational Taxonomic Unit)*:** É uma definição operacional utilizada para classificar grupos de indivíduos próximos. Tipicamente são sequências agrupadas de acordo com sua semelhança uns com os outros (geralmente 97% de similaridade). Se ainda não consegue entender, tente pensar que uma OTU é como uma caixa, pode ter nome (espécie conhecida) ou não (espécie desconhecida) e dentro dela se encontram todos os indivíduos similares (normalmente mínimo 97% de similaridade).
**ASV (*Amplicon Sequence Vatiant*):** É o termo usado para se referir a sequências de DNA individuais. Aqui já não pode pensar que são caixinhas, porque as ASVs não se agrupam por similaridade. Este conceito considera que uma sequência é diferente de outra só com a variação de uma base nucleotídica.

**Mas por qué tem dois conceitos diferentes?** Bom, porque os pesquisadores observaram que nem sempre duas espécies são diferentes porque suas sequências são 3% dissimalares. Tem casos onde tem sequências que são menos de 3% diferentes entre se e são espécies diferentes.

Dito isto, agora vamos a entender a diferença no tipo de clusterização para cada um (ver figura).

 Figura baseada na figura S1 in *Callahan et al. 2016*

A figura mostra uma visão geral do processo de metabarcoding, com os vieses principais que afetam potencialmente a acurácia do sequenciamento. Na amostra **(A)** estão presentes várias espécies com biomassa diferente (indicada pelo tamanho da bolinha) e haplótipos distintos (indicados pela cor). Após a extração de DNA, o gene *barcode* é amplificado usando PCR (B), que pode distorcer a abundância das sequências, e também poderia não amplificas os taxa devido ao viés do iniciador ou à profundidade de sequenciação insuficiente no caso de taxa sub-representada/rara (rosa - sp.5). No processo de high-throughput sequencing (C), muitas novas variantes da sequência falsa são geradas devido a erros de sequenciamento, formação de quimera e mistura de amostras multiplexadas. O impacto desses erros geralmente pode ser reduzido pela rigorosa filtragem de qualidade e pelo agrupamento de seqüências semelhantes em OTU. Normalmente, apenas a sequência mais abundante em uma OTU é considerada e usada para identificar as respectivas espécies (centroides), o que, por sua vez, significa que as informações sobre diversidade genética são perdida. (D) No qiime2 se usam estratégias de ***denoising*** para remover seqüências afetadas por erro de um conjunto de dados e reativar as seqüências de haplótipos reais presentes em uma amostra. O denoising consegue discriminar entre um erro do sequenciamento ou de amplificação na PCR e uma espécie diferente (diferença de no minimo uma base nucleotica). O baseamento é se a sequencia estiver muito representada (abundante) o mais provável é que seja uma sequencia certa, mas se a sequencia estiver em baixa abundância provavelmente é devido a um erro de PCR ou de sequenciamento. Em resumo usando o conceito de OTU, neste exemplo se considerariam só 3 OTUs ou espécies diferentes, enquanto que usando a abordagem de denoising é possivel obter 5 ASVs ou espécies diferentes.

---

**Continuemos com os comandos...**

## 5 Denoising, joining reads e remoção de chimeras com DADA2
DADA2 é um programa que é capaz fazer o *denoising*, unir as sequencias (R1+R2) e eliminar as chimeras. Além, com ele você também poderia fazer o controle de qualidade. Se você quiser explorar todas possibilidades que o DADA2 oferece digite `qiime dada2 --help` 🆘 Dependendo do tipo de sequencias que você tenha deve escolher entre `denoise-paired` (sequencias paired-end illumina), `denoised-pyro` (sequencias single-end por pirosequenciamento) e `denoise-single` (sequencia single-end illumina).
No nosso caso, nossas sequencias são paired-end obtidas em plataforma Illumina.
Você também pode ver todos os parâmetros que oferece o comando `denoised-paired`. Digitando `qiime dada2 denoise-paired --help`.

Por enquanto para nosso exemplo e como a gente já fez o controle de qualidade, não vamos usar os parâmetros para trimagem das sequencias.

## Rodar DADA2
qiime dada2 denoise-paired --i-demultiplexed-seqs 04.ImportedReads/reads_trimmed.qza \
                           --p-trunc-len-f 0 \
                           --p-trunc-len-r 0 \
                           --output-dir 05.Dada2Output

Os arquivos de saída desse comando são
* **denoising_stats.qza**: Este arquivo contém as estatísticas do processo, ou seja, quantas sequencias ficaram depois de unir o R1+R2, remover as chimeras e fazer o *denoising*. Mas ele 🚫**NÃO PODE SER ABERTO**🚫, por ser de tipo `.qza`. Tem que ser transformado a `.tsv`.
* **representative_sequences.qza**: Neste arquivo estão cada uma das sequencias representativas encontradas nas amostras. Para aceder as sequencias é preciso transformar a `.fasta`
* **table.qza**: nesta tabela se encontram cada um dos ASVs e as frequências deles.

Para convertir o arquivo de saída das estatisticas do denoising de `.qza` para `.tsv`

qiime tools export --input-path 05.Dada2Output/denoising_stats.qza --output-path 05.Dada2Output
## Ver o arquivo
nano 05.Dada2Output/stats.tsv
## Completar seu arquivo QC_controle.xlsx!!

Para convertir o arquivo de saída das sequencias representativas`.qza` para `.fasta`

qiime tools export --input-path 05.Dada2Output/representative_sequences.qza --output-path 05.Dada2Output/

O arquivo de saída vai ser chamado de `dna-sequences.fasta`

## 6 Asignar taxonomia aos ASVs
**Lembre-se**:exclamation: Não só se limite a cópiar e colar o comando 🙄, tente entender que está acontecendo e quais são as partes do comando 🤓.

Para asignar a taxonomia tem que usar uma base dados **(db)** "treinada" (adaptada para Qiime2). Nos já temos várias **db** prontas.
Para 👀 quais **db** temos disponivéis:

ls /home/bioinfo/Documentos/Databases/16S/Silva

16S full length (```classifier_silva_132_99_16S.qza```)
16S regiões V3/V4 (`classifier_silva_132_99_16S_V3.V4_341F_805R.qza`)
16S regiões V4/V5 (`classifier_silva_132_99_16S_V4.V5_515F_926R.qza`)
16S regiões V6/V8 (Arqueia) (`classifier_silva_132_99_16S_V6.V8_A956F_A1401R.qza`)
16S regiões V6/V8 (`classifier_silva_132_99_16S_V6.V8_B969F_BA1406R.qza`)
16S regiões V3/V4 (release 138) (`classifier_silva_138_99_16S_V3.V4_341F_805R.qza`)
16S região V4 (release 138) (`classifier_silva_138_99_16S_V4_515F_806R.qza`)

Você vai ver o listado de **db** disponivéis no nosso servidor. Principalmente usamos a **db** Silva que é mais atualizada e acurada. Vai encontrar o *release 132* (04/18) e o *release 138* (12/19) e para diferentes regiões do 16S. **Cuide de escolher bem sua db.**

Agora vamos para o comando de assignação da taxonomia.

## Criar um novo diretório para a classificação taxonômica
mkdir 06.TaxonomyClassification
## Run
qiime feature-classifier classify-sklearn \
  --i-classifier /home/bioinfo/Documentos/Databases/16S/Silva/classifier_silva_138_99_16S_V4_515F_806R.qza \
  --i-reads 05.Dada2Output/representative_sequences.qza \
  --o-classification 06.TaxonomyClassification/taxonomyclassification.qza

Se der este `[Errno 28] No space left on device`. A explicação do erro é que qiime2 enquanto vai rodando comandos vai criando arquivos temporais e a 📁 onde se estão armazenando está lotado, por tanto vamos a criar uma 🆕 📁 para esses temporais (**ISTO SÓ PRECISA SER FEITO SE SAIR E SÓ SERÁ A PRIMEIRA VEZ QUE VOCE FOR RODAR ESTE COMANDO**). Você pode arrumar com os seguintes comandos e instruções:

## Primeiro precisa desativar o ambiente virtual do qiime2
conda deactivate
## crie uma pasta chamada tmp/ CUIDADO crie ela na sua pasta principal /data/usuário
cd ..
mkdir tmp
## Agora vamos a avisar pro qiime qual é a nova pasta pros arquivos temporais. Lembrando que sempre que está a palavra usuário é para ser substituida como o seu nome de usuário
export TMPDIR=/data/usuário/tmp/
## Conferir
echo $TMPDIR
## Tem que mostrar o caminho para a nova pasta criada.
## Agora roda de novo o comando de asignação da taxonomia. 

Converta o arquivo `taxonomyclassification.qza` para `taxonomyclassification.qzv`

qiime tools export \
   --input-path 06.TaxonomyClassification/taxonomyclassification.qza --output-path 06.TaxonomyClassification

## 7 Avaliação com *BLAST*

O desempenho da classificação taxonômica é difícil de avaliar sem uma referência, no entanto é bom você fazer uma verificação básica comparando as atribuições taxonômicas com os principais acertos do BLASTn para determinados ASVs.

Primeiro é necessário transformar o arquivo `representative_sequences.qza` de saída do **Dada2** e abrir no 🔗[Qiime2 View](https://view.qiime2.org).

qiime feature-table tabulate-seqs --i-data 05.Dada2Output/representative_sequences.qza \
                                  --o-visualization 05.Dada2Output/representative_sequences.qzv

representative_sequences.qzv
Este arquivo facilita fazer o "*blast*" de determinados ASVs no banco de dados NCBI nt. Ao comparar esses hits do BLAST com a atribuição taxonômica dos ASVs (etapa 10), você pode ter certeza de que as atribuições taxonômicas funcionaram corretamente. É uma boa idéia selecionar ~ 5 ASVs para BLAST para essa validação, que deve ser de grupos taxonomicamente diferentes, como diferentes filos, de acordo com o classificador taxonômico.

## 8 Eliminar ASVs contaminantes

Depois de asignar a taxonomia aos ASVs, devemos remover os ASVs que são provavelmente contaminantes. Os contaminantes mais comuns nos dados de sequenciamento 16S são as sequências mitocondriais e cloroplásticas, que podem ser removidas excluindo qualquer ASV que contenha esses termos em seu rótulo taxonômico. Você também poderia excluir qualquer ASV que não seja classificado no nível do filo, uma vez que é mais provável que essas sequências sejam quimeras. Ou você também pode **não** eliminar esta informacão, se você estiver estudando um ambiente pouco caracterizado, onde tenha uma boa chance de identificar novos filos.

O comando usado para eliminar sequencias mitocondriais e de cloroplasto é o seguinte:

qiime taxa filter-table \
   --i-table 05.Dada2Output/table.qza \
   --i-taxonomy 06.TaxonomyClassification/taxonomyclassification.qza \
   --p-exclude mitochondria,chloroplast \
   --o-filtered-table 05.Dada2Output/table_filt_contam.qza

Transfor o arquivo `table_filt_contam.qza` para `table_filt_contam_summary.qzv`, e use o Qiime2 View para abrir o arquivo e explorar diferentes informações e comparar com a tabela sem filtrar `table_summary.qzv`

qiime feature-table summarize \
   --i-table 05.Dada2Output/table_filt_contam.qza \
   --o-visualization 05.Dada2Output/table_filt_contam_summary.qzv

Se você comparar o número de ***features*** de cada amostra (aba *Interactive sample detail*) entre a tabela inicial e a filtrada, pode reparar que não tem diferença, por tanto neste caso não tinha contaminação com DNA mitocondrial e/ou de cloroplastos.
Enquanto a filtrar ASVs sem afiliação taxonômica, é melhor primeiro ter uma visão mais detalhada dos resultados obtidos e se for o caso posteriormente eliminar esses ASVs dentro ambiente do **software R**.

## 9 Construir a 🌲 com os *plugins* *FastTree* e *MAFFT*

## Crie um diretório para guardar a árvore
mkdir 07.PhylogeneticTree
## Primeiro você vai fazer um alinhamento múltiplo com MAFFT
qiime alignment mafft --i-sequences 05.Dada2Output/representative_sequences.qza \
                        --p-n-threads 4 \
                        --o-alignment 07.PhylogeneticTree/rep_seqs_aligned.qza
## 
qiime alignment mask --i-alignment  07.PhylogeneticTree/rep_seqs_aligned.qza \
                       --o-masked-alignment 07.PhylogeneticTree/rep_seqs_aligned_masked.qza
         
## Constuir a árvore com FastTree
qiime phylogeny fasttree --i-alignment 07.PhylogeneticTree/rep_seqs_aligned_masked.qza \
                         --p-n-threads 4 \
                         --o-tree 07.PhylogeneticTree/rep_seqs_aligned_masked_tree.qza

## Enraiçar a árvore

qiime phylogeny midpoint-root --i-tree  07.PhylogeneticTree/rep_seqs_aligned_masked_tree.qza \
                              --o-rooted-tree 07.PhylogeneticTree/rep_seqs_aligned_masked_tree_rooted.qza
## Exportar a árvore
qiime tools export /
--input-path 07.PhylogeneticTree/rep_seqs_aligned_masked_tree_rooted.qza \
--output-path 07.PhylogeneticTree/

Agora você criou o arquivo da 🌲 `tree.nwk`. Vamos usa-o em análises mais para frente, onde vamos a calcular as distâncias filogenéticas.

## 10 Metadata
Crie um arquivo de metadata chamado `sample-metadata.txt`. Nesse arquivo você vai colocar o máximo de informação que permita descrever e diferenciar suas amostras.

## Crie o arquivo com o editor de texto nano
nano sample-metadata.txt
## Coloque quantas colunas (variavéis) você quiser que descrevam suas amostras. Aqui o exemplo para o dataset exemplo do tutorial
sample-id	SampleName	StageOfTreatment	Local	Season	
#q2:types	categorical	categorical	categorical	categorical
amostra1	AFA	AnaerobicFilter	Factory	Autumn
amostra2	STW	SepticTank	School	Winter
amostra3	AFS	AnaerobicFilter	Factory	Summer
amostra4	STS	SepticTank	School	Spring

## A segunda fila é para classificar as variavéis segundo sua natureza (p.e. categorical ou numerical) 

Se deve ver assim:

## 11 Barplot 📊
Você pode construir vários gráficos dentro do qiime2, entre eles 📊. Para dar uma primeira olhada na composição taxonômica de seus dados. No entanto para publicações e trabalhos é melhor construir os gráficos com outras ferramentas como o **software R**.

## Crie um diretório para guardar os gráficos
mkdir 08.Graphs
## Make a barplot
qiime taxa barplot \
  --i-table 05.Dada2Output/table.qza \
  --i-taxonomy 06.TaxonomyClassification/taxonomyclassification.qza \
  --m-metadata-file sample-metadata.txt \
  --o-visualization 08.Graphs/taxa_barplots.qzv

📊 on Qiime2 View


## 12 Rarefaction curves 📈
Para construir as curvas de rarefação é importante normalizar as amostras todas ao mesmo tamanho, para isto você pode transformar o arquivo `table.qza` para `table.qzv`, e usando o Qiime2 View pode abrir o arquivo e explorar diferentes informações

qiime feature-table summarize \
           --i-table 05.Dada2Output/table.qza \
           --o-visualization 05.Dada2Output/table_summary.qzv \
           --m-sample-metadata-file sample-metadata.txt


Com este arquivo pode ser visualizado o número de ***features*** (ASVs) finais em cada amostra (aba *interactive sample detail*).

Para este caso a amostra com menor número de *features* é amostra2 com 57785. Então vamos a cortar a essa profundidade para as curvas de rarefação. Com este comando você pode construir as curvas de rarefação usando diferentes mêtricas de alfa diversidade, como "obseved_otus", "chao1", "shannon" e "simpson".

## Crie um diretório para armazenar as curvas de rarefação
mkdir 09.RarefactionCurves
## Rodar
qiime diversity alpha-rarefaction --i-table 05.Dada2Output/table.qza \
                                    --p-max-depth 57785 \
                                    --p-steps 10 \
                                    --p-metrics shannon \
                                    --p-metrics observed_otus \
                                    --p-metrics simpson \
                                    --p-metrics chao1 \
                                    --m-metadata-file sample-metadata.txt \
                                    --o-visualization 09.RarefactionCurves/rarefaction_curves.qzv

Como você pode reparar o arquivo de saída já é de tipo `.qzv` (visualização) pra 👀 no Qiime2 View
rarefaction_curves.qzv.

Todas as amostras atingiram o *plateau* por tanto o esforço amostral foi suficiente para acessar a toda a diversidade microbiana das amostras.

## 13 Análise de alfa e beta diversidade
**PARA!** Antes de continuar é **muito importante** você entender que é alfa e beta diversidade, as diferentes mêtricas usadas para a análises e comparações de microbiomas e até testes estatísdicos adequados para este tipo de dados.

Aqui só alguns links e referências recomendados para você revisar antes de continuar 🤓 🔗[Alpha and beta diversity metrics](http://www.evolution.unibas.ch/walser/bacteria_community_analysis/2015-02-10_MBM_tutorial_combined.pdf)
📖 Magurran AE. Measuring biological diversity: John Wiley & Sons, 2013.
🔗[GUide to STatistical Analysis in Microbial Ecology (GUSTA ME)!](https://sites.google.com/site/mb3gustame/home)

**Aqui um resumão do maisss básico** 🤓

**Alfa diversidade** (Dentro das amostras)

***Riqueza de espécies***, quantas espécies ou ASVs diferentes podemos detectar em cada amostra?"
*Mêtricas:* Espécies observadas (observed_otus), ACE, Chao1.

***Diversidade de espécies*** (Shannon index) "How different?"
Como estão distribuidos (equilibrados) estão os microrganismos entre si? Temos uniformidade de espécies (nível de abundância semelhante) ou algumas espécies dominam mais que outras?
*Mêtricas:* Shannon, Simpson.

**Beta diversity** (Entre as amostra)

Quão diferente é a composição microbiana em um ambiente (amostra) em comparação com outro(a)?

A beta diversidade mostra a diferença entre comunidades microbianas de diferentes ambientes (amostras). O foco principal está na diferença nos perfis de abundância taxonômica de diferentes amostras.

***Mêtricas:***

Bray–Curtis dissimilarity
- Baseada na abundância ou *read count*
- Diferenças das abundâncias microbianas entre duas amostras (p.e. A nível de espécie)
Valores desde 0 a 1
0 significa que ambas amostras compartilham a mesma espécie exatamente na mesma abundância.
1 significa que ambas amostras têm espécies e abundâncias totalmente diferentes.

Jaccard distance
- Baseada na presença ou ausência de espécies (Não leva em conta a informação de abundância)
- Diferenças na composição microbiana entre duas amostras
0 significa que ambas amostras compartilham exatamente as mesmas espécies.
1 significa que ambas amostras não tem nenhuma espécie em comum

UniFrac
- Distâncias filogenéticas (árvore filogenética)
- Baseada no comprimento da ramificação compartilhada entre duas amostras ou exclusiva para uma ou outra amostra.
unweighted UniFrac: baseado só em distâncias filogenéticas (Não leva em conta a informação de abundância)
weighted UniFrac: os comprimentos dos ramos são ponderados por abundância relativa (inclui informações filogenéticas e de abundância)

---
***Continuando...***

Nesta etapa vamos a gerar gráficos e tabelas para interpretar a alfa e beta diversidade do *dataset* analisado.

A profundidade deve ser a mesma usada na construção das curvas de rarefação baseada no arquivo `table_summary.qzv`

qiime diversity core-metrics-phylogenetic --i-table 05.Dada2Output/table.qza \
                                          --p-sampling-depth 57785 \
                                          --m-metadata-file sample-metadata.txt \
                                          --p-n-jobs 4 \
                                          --i-phylogeny 07.PhylogeneticTree/rep_seqs_aligned_masked_tree_rooted.qza \
                                          --output-dir 10.AlphaBetaDiversity
## Para calcular Chao1
qiime diversity alpha --i-table 05.Dada2Output/table.qza \
                        --p-metric chao1 \
                        --o-alpha-diversity 10.AlphaBetaDiversity/chao1_vector.qza
## para calcular Simpson
qiime diversity alpha --i-table 05.Dada2Output/table.qza \
                        --p-metric simpson \
                        --o-alpha-diversity 10.AlphaBetaDiversity/simpson_vector.qza

São varios os arquivos de saída do anterior comando. São de três tipos.
* **Vector**: São os valores de uma mêtrica de alfa diversidade para todas as amostras, organizados em uma tabela. As mêtricas que se organizam em vectores são: observed_otus, shannon, evenness, simpson, chao1, faith_pd.
* **Matrix**: São matrizes de dissimilaridade de diferenter mêtricas, tais como: bray_curtis, jaccard, unweighted_unifrac, weighted_unifrac.

você pode organizar uma tabela com todos os vectores

# Para organizar uma tabela com todos os indices de alfa diverside e estimadores

qiime metadata tabulate --m-input-file sample-metadata.txt \
                          --m-input-file 10.AlphaBetaDiversity/shannon_vector.qza \
                          --m-input-file 10.AlphaBetaDiversity/observed_otus_vector.qza \
                          --m-input-file 10.AlphaBetaDiversity/simpson_vector.qza \
                          --m-input-file 10.AlphaBetaDiversity/chao1_vector.qza \
                          --o-visualization 10.AlphaBetaDiversity/alfadiversidade_all.qzv

Descarrega o arquivo `alfadiversidade_all.qzv` e vai ver 👀 ele no Qiime2 View.


## 14 Exportação dos arquivos
Os arquivos finaís que você irá precisar são:
* `table.qza`, o qual irá se convertir em `Otu_table.tsv`. Por default é chamada assim mas no caso do Qiime2 são ASVs. Esta tabela contém as frequencias de cada ASVs em cada amostra
* `taxonomyclassification.qza`, o qual irá se convertir em `Taxonomy_table.tsv`. Esta tabela contém a taxonomia de cada ASV.
* `tree.nwk`
Os dois primeiros devem ser transformados a `.tsv` para ser usados nas análises seguintes.

## Crie um diretório para salvar os arquivos finais 
mkdir 11.ExportFiles
## Primeiro transformar .qza para .biom
qiime tools export --input-path 05.Dada2Output/table.qza \
                     --output-path 11.ExportFiles/
## .biom para .tsv
biom convert -i 11.ExportFiles/feature-table.biom \
               -o 11.ExportFiles/Otu_Table.tsv \
               --to-tsv \
               --table-type "OTU table"

Abra o arquivo gerado `Otu_Table.tsv file` e troque #OTU ID por OTUID.

## Abrir o arquivo
nano 11.ExportFiles/Otu_Table.tsv
## Depois de trocar o #OTU ID por OTUID, sair
[Ctrl + X]
## Deseja salvar as modificações
[s]
## Nome do arquivo (o mesmo)
[Enter]

Agora exporte a tabela de taxonomia

qiime tools export --input-path  06.TaxonomyClassification/taxonomyclassification.qza \
                     --output-path 11.ExportFiles/taxonomy

Abrir o arquivo gerado `taxonomy.tsv`, troque Feature ID por OTUID

nano 11.ExportFiles/taxonomy/taxonomy.tsv
## Depois de trocar o Feature ID por OTUID, sair
[Ctrl + X]
## Deseja salvar as modificações
[s]
## Nome do arquivo (o mesmo)
[Enter]

Estes arquivos serão usados para nosso tutorial 🗒 de 🔗[ análises de microbiomas usando o *software R*.](https://)

# Final `QC_controle.xlsx`

| SampleID | Raw_reads | Post_trim_primers | Perda até aqui| Post_QC_Trim | Perda até aqui | % Perda | Artefato | filtered | Denoised | merged_reads|non-chimeric | Perda até aqui | % perda total |
| -------- | --------- | ------ | ----------------- | --------- | --------- |:---------:| --------- | ---------- | ------ | -------- | ---------- | ------ | ------------ | ---------- | ---------------- | ------------ |
|sample1|153569|148499|5070|97420|56149|36.56|97420|97420|94322|74928|74691|78878|51.36|
|sample2|120080|115727|4353|69894|50186|41.79|68894|68894|67978|57785|57785|62295|51.88|
|sample3|171262|164997|6265|104537|66725|38.96|104537|104537|101551|77947|77763|93499|54.59|
|sample4|113615|110059|3556|69468|44147|38.86|69468|69468|67875|60075|60075|53540|47.12|

---

## FERRAMENTAS ⚒️ IMPORTANTES PARA AS ANÁLISES DOWNSTREAM
🔗[RAW GRAPHS:](https://rawgraphs.io/) Para construir gráficos, você precisa formatar as tabelas no formato que o site exige.

🔗[Microbiome Analyst:](https://www.microbiomeanalyst.ca/) É uma ⚒️ *on-line* para análises de microbiomas. Você precisará as saídas do qiime2 -`Otu_Table.tsv`, `taxonomy.tsv`, `sample_metadata.txt`, `tree.nwk` (optional)-.

## Algumas referências
📄 **Qiime2:** Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet C, Al-Ghalith GA, et al. QIIME 2: Reproducible, interactive, scalable, and extensible microbiome data science. PeerJ Preprints, 2018.

📄 **FASTQC:** Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010.

📄 **Trimmomatic:** Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 2014; 30.

📄 **Cutadapt:** Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. journal 2011; 17: 10-12.

📄 **DADA2:** CCallahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nature methods 2016; 13: 581.

📄 **SILVA:** Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic acids research 2012; 41: D590-D596.

**FIM** :sparkle: