---
# System prepended metadata

title: EasyMicrobiomes

---

# EasyMicrobiomes

Análise de microbiomas com USEARCH/MicrobiomeAnalyst.

{%slideshare RafaelCorreiadaSilva2/fundamentao-de-microbiomas-com-16s-usearchqiime2 %}

Tecnicamente, esse método funciona para qualquer conjunto de dados, em qualquer máquina Linux, até perto de 40 amostras, seguindo o mesmo protocolo trocando apenas alguns parâmetros!

[toc]

Caso seja a análise de outro marcador (ITS), a única coisa que será trocada são os bancos de referência. [Link para o banco de fungos]

## WinSCP Portátil

- [Link](https://cdn.winscp.net/files/WinSCP-5.15.3-Portable.zip?secure=FqYv-JZ-t5Uuz7HNCM-tYA==,1564093656)
- [Dados: Gabriel](https://drive.google.com/a/unesp.br/uc?id=1X1Frt5nVqSrBjbWQ7oe30qaqo_usFsu9)

## 1. Colocando ferramentas e dados no lugar

Primeiro, vamos baixar o conjunto de programas e dados.
Primeiro, vá em seu `HOME`, e crie a pasta "16S_avengers".
Não esqueça de abrir um `screen`!

```shell=
## Volte para o home
cd

## Crie a pasta 16S_avengers
mkdir 16S_avengers

## Entre na pasta
cd 16S_avengers

## Declara a pasta como 'pasta mestre'
MASTER=$(echo $PWD)
```

Feito isso, iremos então fazer o download dos programas necessários.

### 01A. Download dos programas e dados

Vamos preparar um diretório com o esquema "Buffalo" (análises, dados e scripts) para que tudo fique organizadinho.

O USEARCH pode ser baixado gratuitamente do site: www.drive5.com com o preenchimento do e-mail.

- Arquivo 01: [Scripts utilizados](https://puu.sh/DWQyz.zip). Contém 03 programas: (1) o USEARCH, (2) e (3): scripts feitos por mim para converter a taxonomia num formato aceitável para o [Microbiome Analyst](https://www.microbiomeanalyst.ca).
- Arquivo 02: [Bancos de dados utilizados](https://puu.sh/DWQAU.zip), contém o banco de dados RDP Database, uma base de sequências de 16S extremamente curado, isto é, tem-se certeza de que todas as sequências pertencem aos organismos que deram nome a elas.
- Arquivo 03: [Minha coleção de oligos](https://puu.sh/DWR79.fasta), que servirá para buscar nas sequências sua presença, ausência e tamanho.

Os dados brutos já estão disponíveis no servidor. 

Primeiro, crie as pastas:

```shell=
## Cria as pastas necessárias
mkdir -p dados/{brutos,referencias,oligos}
mkdir analises
mkdir scripts

## Download do banco de referências
cd dados/referencias
wget https://puu.sh/DWQAU.zip
unzip DWQAU.zip
cd ..
cd ..

## Cópia dos dados brutos
cd dados/brutos
cp /data/16S_avengers/data.zip ./
unzip data.zip
cd ..
cd ..

## Cópia dos oligos
cd dados/oligos
wget https://puu.sh/DWR79.fasta
mv DWR79.fasta adapters.fasta
cd ..
cd ..

## Download dos scripts
cd scripts
wget https://puu.sh/DWRnz.zip
unzip DWRnz.zip

## Faz com que todos os scripts se tornem executáveis
chmod 777 ./*

## Retorna pra pasta central
cd ..

## Para ver como está:
tree
```

### 02A. Navegando nas pastas

Ao final de todos os downloads e instalação, o seu `tree` deve mostrar que tudo que você precisa está baixado: as referencias (incluindo oligos), os seus dados brutos, e também os scripts necessários:

![](https://i.imgur.com/zYNO8Ua.png)




### 03A. Preparando os diretórios restantes

Cada análise será feita passo-a-passo, e guardaremos cada um desses passos em sua pasta apropriada, para que depois não façamos confusão entre eles. 

Assim, você pode copiar e colar os comandos abaixo, criando as pastas necessárias:

```shell
## Volte na raíz da pasta mestre
cd $MASTER

## Cria todas as outras pastas necessárias
mkdir -p ./analises/E01_A_fastx_info/logs
mkdir -p ./analises/E01_B_fastq_eestats2/logs
mkdir -p ./analises/E01_C_search_oligodb/logs
mkdir -p ./analises/E02_mergepairs/logs
mkdir -p ./analises/E03_trimming/logs
mkdir -p ./analises/E04_orient/logs
mkdir -p ./analises/E05_filter/logs
mkdir -p ./analises/E06_derrep/logs
mkdir -p ./analises/E07_singletons/logs
mkdir -p ./analises/E08_pick_otus_zotus/logs
mkdir -p ./analises/E09_otu_tables/logs
mkdir -p ./analises/E10_sintax_taxonomy
mkdir -p ./analises/E11_distance_tree
mkdir -p ./analises/E12_microbiome-analyst

## Para conferir seu trabalho, faça novamente um tree:
tree
```

![](https://i.imgur.com/svkLjta.png)




## 2. Pré-processamento

Nessa seção, iremos deixar nossas sequências brutas na melhor qualidade possível, além de deixá-las sem repetições e informações indesejadas para que possamos inferir as OTUs verdadeiras.

Primeiro, é interessante ver o que está acontecendo com nossas sequências

### 02A. `fastx_info` Verificação de qualidade

Esse comando permite ver, textualmente, como suas sequências estão em termos de qualidade por base, tamanho, e outras estatísticas menos úteis.

Vamos fazer esse comando para uma sequência, e depois para todas as sequências:

#### Modo básico: Somente o Iron Man

```shell=
## Análise das duas sequências do Homem de Ferro
./scripts/usearch11 \
    -fastx_info ./dados/brutos/IronMan_R1_001.fastq
    
./scripts/usearch11 \
    -fastx_info ./dados/brutos/IronMan_R2_001.fastq
```

#### Para todos de uma vez!

Para não ficar entendiamente rodando um por vez, podemos usar um mecanismo que executa as instruções:

:::info
O comando abaixo é **pseudocódigo**, isto é, instruções de computador escrito em linguagem humana. Ele não deve ser executado, mas é muito útil para entender o que um código em linguagem de máquina está fazendo!
:::

```pseudocódigo=
para cada SEQ em "dados/brutos", faça:
    
    diga "Análise da sequência: VAR[SEQ]"
    
    usearch11 fastx_info VAR[SEQ] \
        guarde os resultados na pasta analises
            
    pule linha
    pule linha
feito
```

Traduzido em Linux, fica da forma:

```shell=
for Seq in $(ls ./dados/brutos/*.fastq); do

    echo "Analise da sequência: ${Seq}"
    
    ./scripts/usearch11 \
        -fastx_info ${Seq} \
        -output ./analises/E01_A_fastx_info/$(basename ${Seq} .fastq).txt
        
    echo "Fim da análise da sequência: ${Seq}"
    echo
    echo

done
```

#### Modelo de resultados

![](https://i.imgur.com/XJkgdxx.png)


Assim, os resultados ficarão guardados na pasta `analises/E01_A_fastx_info`, para posterior consulta.

### 02B: `fastq_eestats2`, consulta de tamanhos de poda

Agora que demos um primeiro visual na sequência, precisamos decidir o quanto das sequências vamos manter, pois podemos podar pontas de baixa qualidade.

Isto é feito com o comando `fastq_eestats2`.

#### Básico

Novamente, pode ser feito sequência por sequência (`Básico`), ou para todas as sequências ao mesmo tempo (`Avançado`):

```shell=
## Roda para o Hulk!
### Sequência R1
./scripts/usearch11 \
    -fastq_eestats2 ./dados/brutos/Hulk_R1_001.fastq \
    -length_cutoffs 150,300,20 \
    -ee_cutoffs 1.0,2.0,3.0
    
### Sequência R2
./scripts/usearch11 \
    -fastq_eestats2 ./dados/brutos/Hulk_R2_001.fastq \
    -length_cutoffs 150,300,20 \
    -ee_cutoffs 1.0,2.0,3.0
```

#### Avançado, guardando resultados

```shell=
for fqfile in ./dados/brutos/*fastq
do

    ## Capture o 'nome base'
    bn=$(basename $fqfile .fastq)
    fq2=$(echo $fqfile | sed 's/R1/R2/')
    
    head -n 100000 ${fqfile} > ./analises/E01_B_fastq_eestats2/tmp_R1.fq
    head -n 100000 ${fq2} > ./analises/E01_B_fastq_eestats2/tmp_R2.fq
    
    ## Rode o comando para cada sequência
    ./scripts/usearch11 \
        -fastq_eestats2 ./analises/E01_B_fastq_eestats2/tmp_R1.fq \
        -output ./analises/E01_B_fastq_eestats2/${bn}.txt \
        -length_cutoffs 150,300,20 \
        -ee_cutoffs 1.0,2.0,3.0,4.0
        
done
```

#### Modelo de resultados

![](https://i.imgur.com/SFcVtz3.png)

### 02C: Em busca dos oligonucleotídeos perdidos

Antes de iniciar os garroteios nas sequências, precisamos encontrar os adaptadores de 16S. 

Eu rodo esse procedimento para apenas um par de sequências, pois se as sequências vieram **do mesmo sequenciamento**, montadas com **a mesma biblioteca** (99% dos casos), o adaptador estará no mesmo lugar em todas as sequências.

Assim, vamos procurar os adaptadores nas sequências do Thanos:

```shell=
## R1 do Thanos
./scripts/usearch11 \
    -search_oligodb ./dados/brutos/Thanos_R1_001.fastq \
    -db ./dados/oligos/adapters.fasta \
    -strand both \
    -userout ./analises/E01_C_search_oligodb/Thanos_R1.txt \
    -userfields query+target+qstrand+diffs+tlo+thi+trowdots
   
## R2 do Thanos
./scripts/usearch11 \
    -search_oligodb ./dados/brutos/Thanos_R2_001.fastq \
    -db ./dados/oligos/adapters.fasta \
    -strand both \
    -userout ./analises/E01_C_search_oligodb/Thanos_R2.txt \
    -userfields query+target+qstrand+diffs+tlo+thi+trowdots
```

#### Resultado R1

O oligonucleotídeo em `R1` sempre deve ser encontrado, e aqui a informação mais interessante é o seu tamanho:

![](https://i.imgur.com/IKIb82s.png)

#### Resultado R2

O oligo em `R2`, por sua vez, nem sempre pode ser encontrado. Em nosso caso, ele não foi encontrado (foram encontrados oligos de `Ion`, o que não tem nada a ver com nossas sequências, então provavelmente são biológicos)

![](https://i.imgur.com/S62K4gw.png)

Assim, antes de prosseguirmos, vamos guardar algumas informações, que são os **tamanhos de poda**:

```shell
## Sempre podar no mínimo 05, pontinha de baixa qualidade
## Da R1, serão podados os 19 (oligos)
poda_R1=19
## Da R2, serão podados somente 05 (ponta de baixa qualidade)
poda_R2=05
```

## 3. Inferindo as sequências verdadeiras, ou 'OTUs'

### 03A. `fastq_mergepairs`: unindo R1 e R2

Une todas as sequências `R1` e `R2`:

```shell=
## Une todas de uma vez!
./scripts/usearch11 \
    -fastq_mergepairs ./dados/brutos/*R1*.fastq \
    -fastqout analises/E02_mergepairs/merged.fq \
    -relabel @ \
    -fastq_maxdiffs 30 \
    -fastq_pctid 70
```

Você verá o seguinte resultado:

![](https://i.imgur.com/bACVUps.png)

Agora, temos só uma sequência para continuar os trabalhos, a `merged.fq`, então não será mais necessário usar **todas** as sequências!



### 03B. `fastx_truncate`: definindo tamanho fixo

Com as sequências unidas, agora a próxima etapa é um tanto iterativa, isto é, você pode ter que pará-la algumas vezes até ter o resultado correto.

Esse passo faz com que todas as sequências fiquem do mesmo tamanho, e também remove as pontas indesejadas, com oligonucleotídeos ou bases de baixa qualidade.

Tente várias vezes com tamanhos diferentes de truncagem:

```shell=
## Tamanho de truncagem
tamanho_trunc=250

./scripts/usearch11 \
	-fastx_truncate analises/E02_mergepairs/merged.fq \
	-trunclen ${tamanho_trunc} \
	-stripleft ${poda_R1} \
	-stripright ${poda_R2} \
	-fastqout ./analises/E03_trimming/trimmed.fq
```

Idealmente, devem sobrar a maioria das sequências, com o maior tamanho possível.

### 03C. `orient`: orientando as leituras no banco

Esse passo é recomendado pelo Robert Edgar, mas particularmente noto poucos ganhos com ele (1 a 2% das sequências).

Ele orienta todas as leituras em relação a um banco-de-dados, para que elas possam ser melhor anotadas:

```shell=
## Orienta as leituras
./scripts/usearch11 \
    -orient ./analises/E03_trimming/trimmed.fq \
    -db ./dados/referencias/rdp.database.udb \
    -fastqout ./analises/E04_orient/oriented.fq
```

![](https://i.imgur.com/oSUJtvL.png)


### 03D. `fastq_filter`: filtrando por qualidade

Uma das menores qualidades para conseguir publicar sua análise é de $Q=20$, ou seja, 1 base errada a cada 100 bases. Mas com frequência pode ser maior, o que pode ser visto com base em tentativa e erro.

Na prática, isso significa que o valor deve ser menor que 1% do tamanho de sequências truncado. No caso, usarei ainda mais estringente, o valor de 3.

```shell=
## Define o número de erros
erros=03

## Roda o comando de filtragem
./scripts/usearch11 \
	-fastq_filter ./analises/E04_orient/oriented.fq \
	-fastq_maxee ${erros} \
	-fastaout ./analises/E05_filter/filtered.fa
```


![](https://i.imgur.com/OKcHIga.png)


### 03E. `fastx_uniques`: desreplicação

O procedimento de desreplicação remove sequências duplicadas, pois não adicionam nenhum valor no procedimento da análise. Como aqui não temos interesse na contagem, não há problema:

```shell=
## Desreplica, ou seja, remove tudo que é duplos
./scripts/usearch11 \
	-fastx_uniques ./analises/E05_filter/filtered.fa \
	-fastaout ./analises/E06_derrep/derrep.fa \
	-sizeout \
	-relabel Uniq
```

![](https://i.imgur.com/0QsXOG6.png)


### 03F. `sortbysize`: remoção de 'solteirões'

Solteirões ou _singletons_ são aquelas sequências que aparecem apenas uma vez. De acordo com Edgar, elas muito provavelmente são erros de sequenciamento, e devem ser eliminados:

```shell=
./scripts/usearch11 \
	-sortbysize ./analises/E06_derrep/derrep.fa \
	-fastaout ./analises/E07_singletons/nosingle.fa \
	-minsize 2
```

### 03G. `cluster_otus`: agrupamento das OTUs

Finalmente, reduzimos ainda mais nossos dados, agrupando todas as OTUs similares, que receberão, por fim, um mesmo nome. Este é um processo com uma base teórica muito pesada, e muito contestado, mas é completamente aceito em publicações com metataxonômica.

Esse algoritmo também irá remover quimeras, isto é, sequências esdrúxulas composta de dois organismos diferentes, que são inúteis para a análise :poop:

```shell=
./scripts/usearch11 \
	-cluster_otus ./analises/E07_singletons/nosingle.fa \
	-otus ./analises/E08_pick_otus_zotus/otus.fa \
	-relabel Otu 
```

![](https://i.imgur.com/Nzxig6C.png)

### (Adendo) Gerando zOTUs 

```shell=
./scripts/usearch11 \
    -unoise3 ./analises/E07_singletons/nosingle.fa \
    -zotus ./analises/E08_pick_otus_zotus/zotus.fa
```

## 4. `otutab`: Construindo a tabela de contagens

Agora que temos nossas OTUs, vamos mapear as leituras para que possamos contar quanto de cada uma dela está presente em nossa análise.

Nesse ponto, muitas das sequências que foram descartadas (seja por baixa qualidade, tamanho) irá mapear de volta nas OTUs, então aquelas informações não são perdidas.

Esta é uma etapa um tanto mais demorada, então talvez seja interessante aqui ir buscar um :coffee: !

O comando é:

```shell=
## Cria a tabela de otus, otutab.tsv
./scripts/usearch11 \
    -otutab ./analises/E04_orient/oriented.fq \
    -id 0.97 \
    -otus ./analises/E08_pick_otus_zotus/otus.fa \
    -otutabout ./analises/E09_otu_tables/otutab.tsv \
    -mapout ./analises/E09_otu_tables/map_out.tsv    
    
## para zOTUS
./scripts/usearch11 \
    -otutab ./analises/E04_orient/oriented.fq \
    -id 1.00 \
    -otus ./analises/E08_pick_otus_zotus/zotus.fa \
    -otutabout ./analises/E09_otu_tables/zotutab.tsv \
    -mapout ./analises/E09_otu_tables/zmap_out.tsv
```
![](https://i.imgur.com/JeU4Mdd.png)


## 5. Anotação com `Sintax`
`: Dando nomes às sequências

Com as sequências de OTUs em mãos, vamos agora dar nomes a elas, usando-se o algoritmo `SINTAX`!

```shell=
./scripts/usearch11 \
    -sintax ./analises/E08_pick_otus_zotus/otus.fa \
    -db ./dados/referencias/rdp.database.udb \
    -tabbedout ./analises/E10_sintax_taxonomy/annot.txt \
    -strand both \
    -sintax_cutoff 0.9
```

Batendo os olhos nos resultados:

![](https://i.imgur.com/pGgZOUc.png)


## 6. Análise exploratória no MicrobiomeAnalyst (MA)

Vamos mover tudo que será usado para a pasta do Microbiome Analyst, para depois poder jogar na plataforma online!

### 06A. Conversão da Taxonomia

Para ser aceito no MicrobiomeAnalyst, a taxonomia deve estar num formato específico, canônico dos métodos legado de microbioma (QIIME_01, Mothur).

Como não existia ainda uma forma de converter entre o formato do `SINTAX` e este formato, eu criei dois scripts, um Shell e um Perl para fazer tais conversões.

- [parseSintax.sh](https://pastebin.com/0D4h0pXn)
- [processSintax.pl](https://pastebin.com/HiisVDk3)

Esses scripts não precisam ser entendidos, basta saber que eles convertem o formato completo do SINTAX para um formato mais simples.

```shell=
### Exporta a variável MASTER para o ambiente
export MASTER=$(pwd)

### Cria uma cópia da anotação original, que será alterada
cp ./analises/E10_sintax_taxonomy/annot.txt \
    ./analises/E10_sintax_taxonomy/parse.txt
    
### Converte algumas linhas do arquivo parse
sed -i 's/k:Bacteria/d:Bacteria/g' ./analises/E10_sintax_taxonomy/parse.txt
sed -i 's/"//g' ./analises/E10_sintax_taxonomy/parse.txt

### Roda meu script sobre o arquivo parse,
## Gerando para bacterias uma nova taxonomia
./scripts/parseSintax.sh \
    ./analises/E10_sintax_taxonomy/parse.txt \
    bacteria \
    ./analises/E12_microbiome-analyst/03_taxonomia.txt
    
## No final, irá gravar o arquivo de taxonomia na pasta do MicrobiomeAnalyst.
```

### 06B. Cópia da tabela de OTUs

Aqui nada de estranho acontece, ele só copia os dados gerados no [passo 4](https://hackmd.io/XH-JDMqHRiGOB-fuILvWYA?both#4-Construindo-a-tabela-de-contagens) para a pasta do microbiome analyst:

```shell=
sed 's/OTU ID/NAME/g' ./analises/E09_otu_tables/otutab.tsv > \
	./analises/E12_microbiome-analyst/01_otu_table.txt
```

#### Para zOTUs

```shell=
sed 's/OTU ID/NAME/g' ./analises/E09_otu_tables/zotutab.tsv > \
	./analises/E12_microbiome-analyst/01_zotu_table.txt
```


### 06C. Árvore de distância entre amostras

Calcula a distância entre as amostras, muito similar a UPGMA, utilizado para cálculos que necessitam da distância entre OTUs:

```shell=
## I. Distância entre OTUs
./scripts/usearch11 \
	-calc_distmx ./analises/E08_pick_otus_zotus/otus.fa \
	-tabbedout ./analises/E11_distance_tree/otus_mx.txt \
	-maxdist 0.2 \
	-termdist 0.3
    
## II. Cria a árvore newick
./scripts/usearch11 \
	-cluster_aggd ./analises/E11_distance_tree/otus_mx.txt \
    -treeout ./analises/E12_microbiome-analyst/04_tree.nwk \
    -clusterout ./analises/E11_distance_tree/clusters.txt \
    -id 0.80 \
    -linkage min
```

### 06D. Forjando os metadados

Se tudo der certo, a pasta do MicrobiomeAnalyst deve ter os arquivos:

![](https://i.imgur.com/fUMdWnD.png)

Resta agora, o arquivo que associa cada amostra a seus tratamentos, chamada de **Metadados**. Isto pode ser feito no Excel (mais fácil), ou qualquer editor de texto, como bloco de notas ou Vim (mais difícil). No caso, aqui, faremos no Excel/Calc por estar em todos os computadores do curso:

![](https://i.imgur.com/a4lVYu0.png)

:::warning
Salve esse arquivo no formato texto, e não aquele padrão do Excel (xlsx).
O MicrobiomeAnalyst só aceita os formatos `CSV` e `TXT`.
:::

Idealmente, salve na mesma pasta dos outros arquivos `TXT`:

![](https://i.imgur.com/CU3PMJK.png)

O produto final da nossa análise prévia ficará da forma seguinte:

![](https://i.imgur.com/9tm1Spx.png)

### 06E. Parametrização e normalizações

Acesse o site do MicrobiomeAnalyst, www.microbiomeanalyst.ca.
Feito isso, clique na nossa opção de "Microbiomas":

![](https://i.imgur.com/0yHWu9z.png)

A seguir, carregue seus dados conforme abaixo:

![](https://i.imgur.com/OE3rGAU.png)

Nas duas imagens seguir, podemos ver os tamanhos de bibliotecas, além de detalhes numéricos da nossa análise, que frequentemente são copiados em relatórios/publicações:

![](https://i.imgur.com/aRyNGTc.png)

![](https://i.imgur.com/TD0NB2h.png)

Com isso, damos procedimento a análise. Aqui, ele fornece a opção de cortar sequências, seja por abundância ou por baixa variância entre amostras. Como já cortamos da forma que queríamos, vamos deixar tudo em zero para que sejam incluídos todos os dados.

![](https://i.imgur.com/RERooAQ.png)

Nesta última etapa antes da 'análise divertida', são fornecidas opções de normalizações e transformações.

Estas nada mais são do que aplicar uma conversão matemática nas tabelas, e deve variar para alguns métodos.

Para análise exploratória inicial, é interessante que nenhuma transformação seja feita, e os dados sejam somente escalonados em relação ao total, o que não perturba os dados de forma alguma.

![](https://i.imgur.com/YGJ8vdL.png)

Uma das métricas mais aceitas é rarefação, que ainda é aceita em muitas revistas de ecologia microbiana e microbiologia, mas possui muitos problemas e conflitos, como várias outras áreas da microbiômica. 

### 06F. Diversidade alfa

Métrica que, em geral, discute a riqueza de cada ambiente. Pode ser feito testes de médias, mas em geral é usado apenas como diagnostico do ambiente em isolado.

No MicrobiomeAnalyst, ele é compardo com um teste de médias:

![](https://i.imgur.com/7d2sIc6.png)

### 06G. Diversidade beta

Esta é a métrica que permite a comparação entre ambientes, e quando associado a um PCoA ou NMDS, permite o 'mapeamento' das amostras, dizendo o quão similares ou dessimilares eles são:

![](https://i.imgur.com/bkgQseV.png)


### 06H. Testes de média

Teste de média típico (Kruskal-Wallis), que permite fazer afirmativas do tipo "mais abundante/menos abundante"

![](https://i.imgur.com/KHRItSS.png)

Figura do MetacodeR, que tenta exibir de uma só vez todas as diferenças em táxons (figura pouco informativa, mas 'gracinha', deve ser editada para utilizá-la melhor)

![](https://i.imgur.com/1bXBWdB.png)


### 06I. Agrupamento de amostras

Faz o agrupamento em diversas métricas e mapa de calor, exibindo como se distribuem as taxonomias por organismo:

![](https://i.imgur.com/52TwI0c.png)


### 06J. Predição de fatores importantes

O MicrobiomeAnalyst também faz uma análise de discriminante linear, que elege OTUs ou taxa que são muito importantes para caracterização daquele descritor nos metadados.

![](https://i.imgur.com/MNJpbv7.png)


### 06K. Distribuições globais

Figura mais canônica de análise de microbiomas, possui algumas das análises mais típicas de QIIME2, muito publicadas mas que eu não vejo muita finalidade de engenharia:

![](https://i.imgur.com/xePJSqs.png)

## 07. Predição de metagenomas: PICRUSt2

A nova versão do PICRUSt permite a predição funcional a partir dos marcadores taxonômicos, como se fosse uma 'consulta do RG' para descobrir 'o emprego da bactéria'.

Em termos técnicos, o `PICRUSt2` tenta inferir a função das bactérias a partir dos dados:

1. As sequências das bactérias presentes e;
2. Sua contagem (quanto há de cada uma delas).

Aqui, o termo importante é **tenta inferir**, ou seja, tenta prever. Esse é um método com importantes limitações ([leia aqui sobre elas](https://github.com/picrust/picrust2/wiki/Key-Limitations)), mas que consegue com frequência acertar predições e dar algumas evidências das funcionalidades em um dado ambiente.

Para rodar o PICRUST2, é necessário instalá-lo, o que já está feito na máquina. Caso queira instalar na sua, você deve possuir um ambiente `conda`, e rodar o seguinte código:

```shell=
## Já instalado no servidor!

conda create \
    -n picrust2 \
    -c bioconda \
    -c conda-forge picrust2
```

:::	warning
**Essa etapa é bem demorada, então o código fica disponível, mas os resultados já estão prontos. Caso queira 'ver rodar', aviso que não irá terminar em tempo de aula.**
:::

```shell=
## 1. Ligue o ambiente do picrust2
conda activate picrust2

## 2. Volte para a pasta 'master'
cd ${MASTER}

## 3. Execute a pipeline completa do picrust2
picrust2_pipeline.py \
    -s ./E08_pick_otus_zotus/otus.fa\ 
    -i ./E12_microbiome-analyst/01_otu_table.txt \
    -o ./analises/E13_picrust2 \
    -p 5
```

Com essa etapa rodada, navegue para a pasta do PICRUSt, e copie os arquivos de saída para o seu computador:

```shell
## Descomprime os arquivos
gunzip ./analises/picrust2_test/EC_metagenome_out/*
gunzip ./analises/picrust2_test/KO_metagenome_out

## Copie para seu computador os arquivos:
./analises/picrust2_test/EC_metagenome_out/pred_metagenome_unstrat.tsv
./analises/picrust2_test/KO_metagenome_out/pred_metagenome_unstrat.tsv
```
### 07A. PICRUSt2 no MicrobiomeAnalyst

Com as duas planilhas copiadas, uma única alteração precisa ser feita nos arquivos: a linha de cabeçalho principal deve ser trocada para #NAME, em qualquer editor de texto:

#### Antes

![](https://i.imgur.com/NkJ5Syt.png)

#### Depois

![](https://i.imgur.com/xNyqewK.png)

Feito isso, vá ao site do MicrobiomeAnalyst (www.microbiomeanalyst.ca) e acesse a seguinte opção, `Shotgun Data Profiling`:

![](https://i.imgur.com/FPGXaPH.png)

Agora, abra a aba `Upload a gente abundance table`:

![](https://i.imgur.com/UcCa1uI.png)

Preencha, assim, com os dados de KEGG ou EC number, e também os metadados:

![](https://i.imgur.com/1bp4kDZ.png)


---
# Conceitos importantes

Alguns conceitos interessantes de se revisar e manter em mente sempre que for rodar uma análise de microbioma.

## Alfa-diversidade

Na ecologia, a diversidade alfa é atribuida como a diversidade média de espécies em locais ou habitat em escala local. O termo Alfa-diversidade foi introduzido por Robert Whittaker.

Distinguindo os diferentes taxons (otu) presentes em cada amostra. Sendo então a "riqueza de espécies" o número diferente de espécies em cada amostra.

![](https://drive.google.com/a/unesp.br/uc?id=1t3EeD1asFrpAfEYB7ocImxApvLb7BfX-)

## Banco de dados: SILVA, Greengenes, RDP

Banco de dados Silva rRNA, um recurso online para verificar e alinhar sequencias de RNA ribossomal.
O Silva fornece conjunto de dados abrangentes de qualidade e regularmente atualizados em sequências alinhadas de RNA ribossômico pequeno (16S / 18S, SSU) e subunidade grandes (23S / 28S, LSU) para todos os três domínios da vida (Bacteria, Archaea e Eukarya).

![](https://i.imgur.com/aUyJfpl.png)

O Banco de dados Greengenes é dedicada a Bacteria e Archaea. A classificação é baseada na construção automática da árvore de novo e no mapeamento de classificação de outras fontes de taxonomia (principalmente NCBI). A árvore filogenética é construída a partir de seqüências 16S rRNA que foram obtidas de bancos de dados públicos e passaram por uma filtragem de qualidade.

![](https://i.imgur.com/TMTZadK.png)

O Bnco de dados RDP (Ribossomal Database Project), fornece seqüências de rRNA Bacterial e Archaeal 16S, com controle de qualidade, alinhadas e anotadas, e sequências de rRNA Fúngicas 28S e um conjunto de ferramentas de análise para a comunidade científica.
nformações sobre classificação taxonômica para Bacteria e Archaea são baseadas nos roteiros taxonômicos de Bergey’s manual e LPSN. nformações taxonômicas para fungos são obtidas a partir de uma classificação feita à mão dedicada à taxonomia de fungos.

![](https://i.imgur.com/s4FiIzZ.png)

## Beta-diversidade

Na ecologia Beta diversidade é caracterizada como a diversidade microbiana da comunidade entre duas amostras (dois ambientes). A beta diversidade é a razão entre a diversidade de espécies regionais e locais.

A beta diversidade também foi introduzida por Robert Whittaber.

![](https://drive.google.com/a/unesp.br/uc?id=1XcLlNNZoAyntHdKKTYO_FZhbzdHOe4I1)

## Grampo PNA (_PNA Clamp_)

O Grampo PNA é uma tecnologia baseada em uma PCR com um grampeador PNA, onde esse grampo seleciona e amplifica apenas as sequência alvo de DNA mutada como uma porção menor na mistura com as principais seqüências de DNA do tipo selvagem.  Dessa forma impede a amplificação de DNA mitocondrial e plasmidial para amostras de plantas e tecidos.

![](https://drive.google.com/a/unesp.br/uc?id=1mK5ANE7hixBwSHL3rYQWwCh5h4spdRKx)

## OTUs

A Unidade Taxonômica Operacional (OTUs) é uma definição operacional usada para classificar grupos de indivíduos intimamente relacionados.

Nesse sentido a OTU é uma definição para agrupar indivíduos por semelhanças equivalente, mas não necessariamente, à definição de taxonomia de Lineu ou de toxonomia evolutiva moderna.

![](https://drive.google.com/a/unesp.br/uc?id=14dLD0joCFpwW6gP-KJxaOVgYkyZhANZK)

Em termos práticos (nos programas), uma OTU é uma sequência que corresponde a um organismo único na análise.

## Pré-processamento

Tudo que acontece antes da inferência de informações biológicas, etapas necessárias antes de recuperar informações.

Os diversos sequenciamentos implicam em erros e deficiências nas sequências-alvo, e geralmente se refere a pré-processamento as etapas que removem ou diminuem tais erros, deixando as sequências adequadas para continuidade da análise.

## Qualidade Phred

O algoritmo Phred é o escore de qualidade sendo uma medida de qualidade e de identificação das bases geradas pelo sequenciamento automático de DNA. As pontuações de qualidade do Phred são atribuídas a cada chamada de base de nucleotídeos.

Ele identifica as extremidade dos dados brutos estabelece o padrão de cada um dos quatro sinais de fluorescência, após isso ele utiliza o método do processamento dos sinais para suavizar as estimativas do dado informacional, por fim é gerado uma seqüência de
bases com qualidades associadas que estará presente no arquivo fastq.

![](https://i.imgur.com/ClAIAuq.png)

O Phred score médio de uma sequência é calculado da seguinte forma:

$$
Q = -10*log(P)
$$

Onde 

- $Q$ é a qualidade e
- $P$ é a probabilidade de erro na sequência.

A maioria das revistas aceita uma qualidade mínima de Q20, que quer dizer no máximo 1 erro a cada 100 bases:

$$
Q = -10*log(0.01) \\
Q = (-10)*(-2) \\
Q = 20
$$

## Rarefação

É uma forma de normalizar experimentos. Num sequenciamento NGS, raramente as amostras possuem a mesma quantidade de sequências, podendo ter disparidades extremas (10 a 20x!).

A rarefação é uma divisão: você divide todas as bibliotecas por um valor que as deixe como se fossem 'iguais'.

Supondo que você tenha as seguintes amostras:

| Amostra | N. de seqs | Fator de rarefação |
| -- | -- | -- |
| A01 | 300 | $A01/A03 = 300/100 = 3$ | 
| A02 | 200 | $A02/A03 = 200/100 = 2$ |
| A03 | 100 | $A03/A03 = 1$ |

Para rarefazê-las, teremos que dividir toda a tabela pela razão entre a amostra e a menor amostra.

Isso é considerável inadmissível para alguns, pois você está intencionalmente jogando sequências fora, sendo que existem métodos estatísticos (no MicrobiomeAnalyst, inclusive) que permite facilmente a comparação entre bibliotecas de diferentes tamanhos.

## Regiões intergênicas eucarióticas (ITS)

As regiões intergênicas eucarióticas (ITS), são espaçadores de trancritos internos do conjunto de genes do RNA ribossomal, denominados de ITS. São marcadores nucleares frequentemente utilizados para análises filogenéticas de diversos grupos eucarióticos incluindo a maioria das famílias das plantas.

![](https://i.imgur.com/SFhZAHp.png)


## Sequência Singleton/"Solteirona"

Um singleton é uma leitura com uma sequência que está presente exatamente uma vez, ou seja, é única entre as leituras.

Se os erros do sequenciadores forem independentes e distribuídos aleatoriamente, é improvável que a seqüência de uma read(leitura) incorreta seja reproduzida por acaso e a maioria dos singletons conterá pelo menos um erro. 

Se é improvável que sequências ruins sejam reproduzidas por acaso, então, se a sequência de uma read(leitura) for encontrada duas ou mais vezes, é provável que esteja correta, embora possa ser a sequência correta de um produto de amplificação ruim, por exemplo, uma quimera.

As leituras que são singletons após a filtragem de qualidade e pre-processamento são, portanto, descartadas e as leituras com abundância de duas ou mais são usadas como entrada para o agrupamento de OTU.

## Tabela de OTUs (OTU Table)

Conta quanto de cada sequência única (OTU) está presente no seu sequenciamento, isto é, "quantas vezes cada organismo ocorre". Sempre tem o seguinte formato:

![](https://i.imgur.com/xDjEe5h.png)


# Links interessantes/úteis

- [Documentação do USEARCH](www.drive5.com)
- [Artigo: Por que rarefazer é inadmissível](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003531)
- [Crítica: abundância de leituras não tem correlação alguma com a abundância real da comunidade](https://drive5.com/usearch/manual/amplification_bias.html)
- [Crítica: os bancos de dados SILVA/Greengenes possuem muitos erros devido às predições computacionais](https://drive5.com/usearch/manual/faq_tax_db.html)
- [O problema com Singletons](https://drive5.com/usearch/manual/singletons.html)
