# Control de qualidade e trimagem
Aqui você vai aprender 🤓 a analisar a qualidade de sequências obtidas por sequenciamento em plataforma Illumina e a *trimar* (filtrar) aquelas ***reads*** (sequências) com baixa qualidade para serem descartadas e melhorar a acurácia dos processos ***downstream***, levando para frente só ***reads*** com alta qualidade.
Para este tutorial 🗒 será usado um conjunto de quatro amostras de sequenciamento do gene 16S rRNA (parte I), quatro amostras de sequenciamento do *Internal transcribed spacer* - ITS (parte II) ou três *sub-samples* de amostras de sequenciamento *shotgun* (parte III)
## 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: +5519981721510
:mailbox_with_mail: Email: khidalgo@javeriana.edu.co
* :male-scientist::male-technologist: Victor Borin Centurion
Biomédico
Doutorando em Genética e Biologia Molecular
Instituto de Biologia - UNICAMP
:iphone:Whastapp: +5519982349780
:mailbox_with_mail: Email: vborincenturion@yahoo.com.br
* :male-scientist::male-technologist: Dr. Tiago Palladino Delforno
Biólogo
:mailbox_with_mail: Email: tiago.palladino@gmail.com
---
## Requisitos
1. Não preparou seu pc 💻 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. Para entender um pouco da estrutura do gene 16S e de barcoding 🔗[vai aqui](https://).
5. Por último 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 ai vai 🔗[aqui](https://) (Obrigatório).
## Alguns tip 💁🏻♀️ interessantes antes de começar
* Lembra da tecla [Tab] :keyboard:? para autocompletar nomes de pastas e/ou arquivos e evitar erros de digitação.
* Este é novo. Com a seta para acima ⬆️ e para abaixo ⬇️ 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.
---
# PARTE I Gene 16S rRNA
## Vamos lá! :beginner:
### Ingressar ao servidor. Instruções 🔗 [aqui](https://)
Estando na sua 📁 `/data/usuário/` vamos criar uma nova 📁 que vai chamar `16S_tutorial`
**Entrada**
```coffeescript=
## Confira onde você está
pwd
```
**Saída**
```coffeescript=
/data/usuário
```
**Entrada**
```coffeescript=
## Crie uma nova pasta
mkdir 16S_tutorial
## confira
ls
```
**Saída**
```coffeescript=
16S_tutorial
```
**Entrada**
```coffeescript=
## Entre na nova pasta
cd 16S_tutorial/
```
**Todos os comandos vão ser rodados estando dentro desta 📁**
Agora vamos a criar uma 📁 para colocar os dados brutos
```coffeescript=
mkdir 00.RawData
```
### Example dataset
O *example dataset* se encontra numa 📁 no servidor. Vamos a criar um 🔗**link** do dataset exemplo na 📁 que você criou. Imagine que é como um atalho dos que se criam no Windows. Assim, evita-se informação repetida e uso de espaço desnecessariamente.
**Entrada**
```coffeescript=
## Cópiar
ln -s /home/bioinfo/Documentos/examples_datasets/16S_v3_v4/* 00.RawData
## Lembre-se que o simbolo * siginifica qualquer caracter. Nesta linha vamos cópiar todos os arquivos dentro da pasta ~/examples_datasets/16S_v3_v4/ para a pasta ~/00.RawData/
## Confira
ls 00.RawData/
```
**Saída**
```coffeescript=
amostra1_1.fq.gz amostra1_2.fq.gz amostra2_1.fq.gz amostra2_2.fq.gz amostra3_1.fq.gz amostra3_2.fq.gz amostra4_1.fq.gz amostra4_2.fq.gz
```
**Nota:** Você poderia trabalhar diretamente com seus dados se você tiver, mas a recomendação 🙂 é primeiro praticar com este *dataset*.
### 1 Trimar ✂️ primers con cutadapt
Para cortar ✂️ os primers das sequências usamos um programinha chamado *cutadapt* instalado dentro do ambiente virtual qiime2-2019.10. (Não sabe o que é um ambiente virtual? Então você não passou pelo tutorial de Conda, 🔗[Vai lá](https://))
*Cutadapt* mostra as *reads* que não começam com primers e removem os primers das *reads* que tem eles.
```coffeescript=
## Active o ambiente virtual
source /opt/Miniconda3/bin/activate qiime2-2019.10
## Crie um novo diretório para armazenar as amostras sem primers
mkdir 01.PrimersTrim
```
A continuação alguns primers padrão usados dependendo da região sequenciada
**Primers para V6-V8**
```
-g ACGCGHNRAACCTTACC \
-G ACGGGCRGTGWGTRCAA \
```
**Primers para V4-V5**
```
-g GTGYCAGCMGCCGCGGTAA \
-G CCGYCAATTYMTTTRAGTTT \
```
**Primers para V3-V4**
```
-g CCTAYGGGRBGCASCAG \
-G GGACTACNNGGGTATCTAAT \
```
Os primers (515F-806R) usados nestas amostras foram, específicos pra V4:
```
-g GTGYCAGCMGCCGCGGTAA \
-G GGACTACNVGGGTWTCTAAT \
```
A continuação o comando, você tera que rodar ele para cada amostra, trocando os nomes das amostras de entrada e dos arquivos de saída
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GTGYCAGCMGCCGCGGTAA \
-G GGACTACNVGGGTWTCTAAT \
-o 01.PrimersTrim/amostra1_r1.fq.gz \
-p 01.PrimersTrim/amostra1_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample1_cutadapt_log.txt' \
::: 00.RawData/amostra1_1.fq.gz ::: 00.RawData/amostra1_2.fq.gz
```
o comando para a amostra 2, note que no flag 🚩 -o e -p os nomes de saída mudaram para `amostra2_r1.fq.gz` e `amostra2_r2.fq.gz` respectivamente. O arquivo 📄 `log.txt` mudou para sample2 e por último os caminhos de entrada mudaram para `00.RawData/amostra2_1.fq.gz` e `00.RawData/amostra2_2.fq.gz`. Assim, você já sabe tudo o que tem que ser mudado cada vez que vai rodar uma amostra diferente. **Preste muita atenção a esto porque você poderia sobre escrever arquivos sem perceber.**:exclamation:
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GTGYCAGCMGCCGCGGTAA \
-G GGACTACNVGGGTWTCTAAT \
-o 01.PrimersTrim/amostra2_r1.fq.gz \
-p 01.PrimersTrim/amostra2_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample2_cutadapt_log.txt' \
::: 00.RawData/amostra2_1.fq.gz ::: 00.RawData/amostra2_2.fq.gz
```
Amostra 3
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GTGYCAGCMGCCGCGGTAA \
-G GGACTACNVGGGTWTCTAAT \
-o 01.PrimersTrim/amostra3_r1.fq.gz \
-p 01.PrimersTrim/amostra3_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample3_cutadapt_log.txt' \
::: 00.RawData/amostra3_1.fq.gz ::: 00.RawData/amostra3_2.fq.gz
```
Amostra 4
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GTGYCAGCMGCCGCGGTAA \
-G GGACTACNVGGGTWTCTAAT \
-o 01.PrimersTrim/amostra4_r1.fq.gz \
-p 01.PrimersTrim/amostra4_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample4_cutadapt_log.txt' \
::: 00.RawData/amostra4_1.fq.gz ::: 00.RawData/amostra4_2.fq.gz
```
Confira que os arquivos :page_facing_up: de saída foram criados
```coffeescript=
ls 01.PrimersTrim/
```
**Saída**
```coffeescript=
amostra1_r1.fq.gz amostra2_r1.fq.gz amostra3_r1.fq.gz amostra4_r1.fq.gz sample1_cutadapt_log.txt sample3_cutadapt_log.txt
amostra1_r2.fq.gz amostra2_r2.fq.gz amostra3_r2.fq.gz amostra4_r2.fq.gz sample2_cutadapt_log.txt sample4_cutadapt_log.txt
```
Agora vamos a usar um pacote de *scripts* que chama microbiome_helper, para criar um *resumo* de todos os arquivos 📄 `log.txt`, para analizar quantas sequências ficaram depois de trimar ✂️ os primers.
```coffeescript=
/home/bioinfo/Documentos/microbiome_helper/parse_cutadapt_logs.py -i 01.PrimersTrim/*log.txt
## Abra o arquivo de texto
nano cutadapt_log.txt
```
**Saída**
```coffeescript=
sample pairs_in read1_match read2_match pairs_passed
sample1 153569 151923 150100 148499
sample2 120080 118814 116955 115727
sample3 171262 169456 166745 164997
sample4 113615 112400 111233 110059
```
Como pode observar, por exemplo da amostra1 entraram 153569 *reads* e sairam 148499.
```coffeescript=
## Remova todos os arquivos .txt intermediarios
rm 01.PrimersTrim/*_log.txt
```
É importante levar um controle das *reads* que vão ficando depois de cada passo.
Para isto crie um arquivo 📄 `.xlsx` (excel), com as seguintes colunas
***Nota***: esta tabela tem colunas que serão usadas no tutorial 🗒 QIIME2 - gene 16S rRNA
| 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 |
| -------- | --------- | ------ | ----------------- | --------- | --------- |:---------:| --------- | ---------- | ------ | -------- | ---------- | ------ | ------------ | ---------- | ---------------- | ------------ |
Você pode encontrar um template desta planilha aqui:
```coffeescript=
/home/bioinfo/Documentos/examples_datasets/QC_controle.xlsx
```
Faça download 🔽 dela no seu computador 💻 usando filezilla. Precisa ajúda para isso? 🔗[Vai aqui](https://)
### 2 Explorar a qualidade das sequências
```coffeescript=
## Crie um novo diretório para armazenar os reportes de qualidade
mkdir 02.FastqcReports
```
Para analisar 🤓 a qualidade das sequências vamos usar um programa chamado Fastqc, ele faz um scan de cada sequência e determina a qualidade delas segundo vários parâmetros, criando um reporte com interfaz gráfica.
```coffeescript=
## Fastqc está instalado no ambiente virtual bioinfo
source /opt/Miniconda3/bin/activate bioinfo
## Rodar fastqc
fastqc -t 10 01.PrimersTrim/* -o 02.FastqcReports/
```
o flag 🚩 -t significa quantos núcleos ou threads você quer que o programa use. Nosso servidor 🖥 tem 22. 10 é bastante para este processo.
Na anterior linha de comando você processou todas as amostras de uma vez só, simbolizado no caracter *. Lea o comando assim: fastqc por favor me mostra a qualidade das sequências de todas as amostras que se encontram dentro da 📁 `01.PrimersTrim/` e coloque os reportes na 📁 `02.FastqcReports/`
Opcionalmente, você pode gerar um reporte único com todas as amostras, usando a ferramenta multiqc
```coffeescript=
multiqc 02.FastqcReports/* -o 02.FastqcReports/
```
Agora para conseguir observar os reportes vai no filezilla e faça download 🔽 deles no seu 💻. Só descarregue os arquivos `.html`
Você vai encontrar dois arquivos por cada amostra o r1 e o r2.
Por enquanto preocupe-se com o parâmetro **Per base sequence quality** que dever ser **Q >30**
Você vai ver esse gráfico de barras :bar_chart:. Lembra no tutorial de tipos de arquivos que tem a explicação do que é *Phred score*? para refrescar a memória, 🔗[vai lá](https://) **é importante**:exclamation:.
O r2 sempre tem uma qualidade mais baixa.
Aqui neste reporte você também pode confirmar o número de sequências (**Basic Statistics**) que foram analisadas, que deve corresponder à saída do cutadapt
Se quiser conhecer mais sobre este reporte, saber que significam os demais parâmetros. 🔗[Vai aqui](https://dnacore.missouri.edu/PDF/FastQC_Manual.pdf).
### 3 Filtragem das *reads* com baixa qualidade
```coffeescript=
## Crie um novo diretório para armazenar as sequências limpas
mkdir 03.CleanData
## Crie outra pasta onde vai armazenar o que NÃO precisa, o lixo
mkdir unpaired
```
Para a trimagem ✂️ das *reads* com baixa qualidade vai usar o programa Trimmomatic
A continuação alguns parâmetros:
**LEADING** Remove as bases com baixa qualidade do começo da sequência. Especifique a qualidade minima para conservar a base.
**TRAILING** Remove as bases com baixa qualidade do final da sequência. Especifique a qualidade minima para conservar a base.
**SLIDINGWINDOW** Cria uma janela deslizante de determinado número de bases. Elimina a janela completa se a média de qualidade é menor da estipulada. WindowSize:requiredQuality.
**MAXINFO** Executa um aro de qualidade adaptável, equilibrando os benefícios de reter as *reads* com o custo de reter bases com erros. targetLength:strictness.
**MINLEN** Remove *reads* com comprimento menor a o estipulado.
Para mais detalhes consulte o 🔗[manual do Trimmomatic](http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf).
```coffeescript=
## Com este comando você vai processar todas as amostras
for i in 01.PrimersTrim/*1.fq.gz
do
BASE=$(basename $i 1.fq.gz)
trimmomatic PE -threads 12 $i 01.PrimersTrim/${BASE}2.fq.gz 03.CleanData/${BASE}1_paired.fq.gz unpaired/${BASE}1_unpaired.fq.gz 03.CleanData/${BASE}2_paired.fq.gz unpaired/${BASE}2_unpaired.fq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MAXINFO:80:0.5 MINLEN:150
done
```
**Atenção!**:exclamation: O parêmetros deste comando podem ser totalmente modificados em função de seus dados e do tipo de sequenciamento (amplicon ou *shotgun*)
O comando deu certo 👍🏼: quando aparece assim:
Você pode reparar que nas quatro amostras sobreviveram aproximadamente 60% das sequências. Você tem que colocar na balança, poucas reads com boa qualidade ou muitas reads com baixa qualidade. Uma boa porcentagem de sequências para ser levado pra frente é minimo de 60% 👍🏼, melhor se fosse mais de 70% 👌🏼.
Mas... vamos a ver 👀 de novo a qualidade como ficou com o fastqc
```coffeescript=
fastqc -t 10 03.CleanData/* -o 02.FastqcReports/
multiqc 02.FastqcReports/*paired* -o 02.FastqcReports/
## Fazer download desde Filezilla
```
Melhorou né? 👍🏼
O r2 também melhorou, no entanto não todas estão na caixa verde, mas é aceptável assim
Olhe na seção **Basic Statitics** o número de sequências restantes para atualizar sua planilha.
Até aqui aprendimos como ver 👀 gráficamente a qualidade das *reads* e a filtrar aquelas com baixa qualidade, para levar as melhores para os análises *downstream*.
# PARTE II ITS
## Vamos lá! :beginner:
### Ingressar ao servidor. Instruções 🔗 [aqui](https://)
Estando na sua 📁 `/data/usuário/` vamos criar uma nova 📁 que vai chamar `its_tutorial`. Lembrando que a palavra usuário vai ser substituida por seu nome de usuário do servidor.
**Entrada**
```coffeescript=
## Confira onde você está
pwd
```
**Saída**
```coffeescript=
/data/usuário
```
**Entrada**
```coffeescript=
## Crie uma nova pasta
mkdir its_tutorial
## confira
ls
```
**Saída**
```coffeescript=
its_tutorial
```
**Entrada**
```coffeescript=
## Entre na nova pasta
cd its_tutorial/
```
**Todos os comandos vão ser rodados estando dentro desta 📁**
Agora vamos a criar uma 📁 para colocar os dados brutos
```coffeescript=
mkdir 00.RawData
```
### Example dataset
O *example dataset* se encontra numa 📁 no servidor. São quatro amostras de sequenciamento do ITS2. Vamos a criar um 🔗**link** do dataset exemplo na 📁 que você criou. Imagine que é como um atalho dos que se criam no Windows. Assim, evita-se informação repetida e uso de espaço desnecessariamente.
**Entrada**
```coffeescript=
## Cópiar
ln -s /home/bioinfo/Documentos/examples_datasets/its_ITS2/* 00.RawData
## Lembre-se que o simbolo * siginifica qualquer caracter. Nesta linha vamos cópiar todos os arquivos dentro da pasta ~/examples_datasets/its_ITS2/ para a pasta ~/00.RawData/
## Confira
ls 00.RawData/
```
**Saída**
```coffeescript=
amostra1_1.fq.gz amostra1_2.fq.gz amostra2_1.fq.gz amostra2_2.fq.gz amostra3_1.fq.gz amostra3_2.fq.gz amostra4_1.fq.gz amostra4_2.fq.gz
```
**Nota:** Você poderia trabalhar diretamente com seus dados se você tiver, mas a recomendação 🙂 é primeiro praticar com este *dataset*.
### 1 Trimar ✂️ primers con cutadapt
Para cortar ✂️ os primers das sequências usamos um programinha chamado *cutadapt* instalado dentro do ambiente virtual qiime2-2019.10. (Não sabe o que é um ambiente virtual? Então você não passou pelo tutorial de Conda, 🔗[Vai lá](https://))
*Cutadapt* mostra as *reads* que não começam com primers e removem os primers das *reads* que tem eles.
```coffeescript=
## Active o ambiente virtual
source /opt/Miniconda3/bin/activate qiime2-2019.10
## Crie um novo diretório para armazenar as amostras sem primers
mkdir 01.PrimersTrim
```
A continuação alguns primers padrão usados dependendo da região sequenciada
**Primers para ITS1**
```
## Primers ITS1-ITS2
-g TCCGTAGGTGAACCTGCGG \
-G GCTGCGTTCTTCATCGATGC \
```
```
## Primers ITS1F-ITS2
-g CTTGGTCATTTAGAGGAAGTAA \
-G GCTGCGTTCTTCATCGATGC \
```
```
## Primers ITS5-ITS2
-g GGAAGTAAAAGTCGTAACAAGG \
-G GCTGCGTTCTTCATCGATGC \
```
**Primers para ITS2**
```
## Primers ITS3-ITS4
-g GCATCGATGAAGAACGCAGC \
-G TCCTCCGCTTATTGATATGC \
```
**Primers para ITS1-5.8S-ITS2**
```
# Primers ITS1-ITS4
-g TCCGTAGGTGAACCTGCGG \
-G TCCTCCGCTTATTGATATGC \
```
Para estas amostras foram usados os primers ITS3-ITS4, específicos para o ITS2:
```
-g GCATCGATGAAGAACGCAGC \
-G TCCTCCGCTTATTGATATGC \
```
A continuação o comando, você tera que rodar ele para cada amostra, trocando os nomes das amostras de entrada e dos arquivos de saída
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GCATCGATGAAGAACGCAGC \
-G TCCTCCGCTTATTGATATGC \
-o 01.PrimersTrim/amostra1_r1.fq.gz \
-p 01.PrimersTrim/amostra1_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample1_cutadapt_log.txt' \
::: 00.RawData/amostra1_1.fq.gz ::: 00.RawData/amostra1_2.fq.gz
```
o comando para a amostra 2, note que no flag 🚩 -o e -p os nomes de saída mudaram para `amostra2_r1.fq.gz` e `amostra2_r2.fq.gz` respectivamente. O arquivo 📄 `log.txt` mudou para sample2 e por último os caminhos de entrada mudaram para `00.RawData/amostra2_1.fq.gz` e `00.RawData/amostra2_2.fq.gz`. Assim, você já sabe tudo o que tem que ser mudado cada vez que vai rodar uma amostra diferente. **Preste muita atenção a esto porque você poderia sobre escrever arquivos sem perceber.**:exclamation:
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GCATCGATGAAGAACGCAGC \
-G TCCTCCGCTTATTGATATGC \
-o 01.PrimersTrim/amostra2_r1.fq.gz \
-p 01.PrimersTrim/amostra2_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample2_cutadapt_log.txt' \
::: 00.RawData/amostra2_1.fq.gz ::: 00.RawData/amostra2_2.fq.gz
```
Amostra 3
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GCATCGATGAAGAACGCAGC \
-G TCCTCCGCTTATTGATATGC \
-o 01.PrimersTrim/amostra3_r1.fq.gz \
-p 01.PrimersTrim/amostra3_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample3_cutadapt_log.txt' \
::: 00.RawData/amostra3_1.fq.gz ::: 00.RawData/amostra3_2.fq.gz
```
Amostra 4
```coffeescript=
parallel --link --jobs 4 \
'cutadapt \
--pair-filter any \
--no-indels \
--discard-untrimmed \
-g GCATCGATGAAGAACGCAGC \
-G TCCTCCGCTTATTGATATGC \
-o 01.PrimersTrim/amostra4_r1.fq.gz \
-p 01.PrimersTrim/amostra4_r2.fq.gz \
{1} {2} \
> 01.PrimersTrim/sample4_cutadapt_log.txt' \
::: 00.RawData/amostra4_1.fq.gz ::: 00.RawData/amostra4_2.fq.gz
```
Confira que os arquivos 📄 de saída foram criados
```coffeescript=
ls 01.PrimersTrim/
```
**Saída**
```coffeescript=
amostra1_r1.fq.gz amostra2_r1.fq.gz amostra3_r1.fq.gz amostra4_r1.fq.gz sample1_cutadapt_log.txt sample3_cutadapt_log.txt
amostra1_r2.fq.gz amostra2_r2.fq.gz amostra3_r2.fq.gz amostra4_r2.fq.gz sample2_cutadapt_log.txt sample4_cutadapt_log.txt
```
Agora vamos a usar um pacote de *scripts* que chama microbiome_helper, para criar um *resumo* de todos os arquivos 📄 `log.txt`, para analizar quantas sequências ficaram depois de trimar ✂️ os primers.
```coffeescript=
/home/bioinfo/Documentos/microbiome_helper/parse_cutadapt_logs.py -i 01.PrimersTrim/*log.txt
## Abra o arquivo de texto
nano cutadapt_log.txt
```
**Saída**
```coffeescript=
sample pairs_in read1_match read2_match pairs_passed
sample1 148237 146743 146351 145189
sample2 140739 139219 138935 137824
sample3 169466 167162 167283 166010
sample4 157229 155430 155029 153877
```
Como pode observar, por exemplo da amostra1 entraram 148237 *reads* e sairam 145189.
```coffeescript=
## Remova todos os arquivos .txt intermediarios
rm 01.PrimersTrim/*_log.txt
```
É importante levar um controle das *reads* que vão ficando depois de cada passo.
Para isto crie um arquivo 📄 `.xlsx` (excel), com as seguintes colunas
***Nota***: esta tabela tem colunas que serão usadas no tutorial 🗒 QIIME2 - gene 16S rRNA
| 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 |
| -------- | --------- | ------ | ----------------- | --------- | --------- |:---------:| --------- | ---------- | ------ | -------- | ---------- | ------ | ------------ | ---------- | ---------------- | ------------ |
Você pode encontrar um template desta planilha aqui:
```coffeescript=
/home/bioinfo/Documentos/examples_datasets/QC_controle.xlsx
```
Faça download 🔽 dela no seu computador 💻 usando filezilla. Precisa ajúda para isso? 🔗[Vai aqui](https://)
### 2 Explorar a qualidade das sequências
```coffeescript=
## Crie um novo diretório para armazenar os reportes de qualidade
mkdir 02.FastqcReports
```
Para analisar 🤓 a qualidade das sequências vamos usar um programa chamado Fastqc, ele faz um scan de cada sequência e determina a qualidade delas segundo vários parâmetros, criando um reporte com interfaz gráfica.
```coffeescript=
## Fastqc está instalado no ambiente virtual bioinfo
source /opt/Miniconda3/bin/activate bioinfo
## Rodar fastqc
fastqc -t 10 01.PrimersTrim/* -o 02.FastqcReports/
```
o flag 🚩 -t significa quantos núcleos ou threads você quer que o programa use. Nosso servidor 🖥 tem 22. 10 é bastante para este processo.
Na anterior linha de comando você processou todas as amostras de uma vez só, simbolizado no caracter *. Lea o comando assim: fastqc por favor me mostra a qualidade das sequências de todas as amostras que se encontram dentro da 📁 `01.PrimersTrim/` e coloque os reportes na 📁 `02.FastqcReports/`
Opcionalmente, você pode gerar um reporte único com todas as amostras, usando a ferramenta multiqc
```coffeescript=
multiqc 02.FastqcReports/* -o 02.FastqcReports/
```
Agora para conseguir observar os reportes vai no filezilla e faça download 🔽 deles no seu 💻. Só descarregue os arquivos `.html`
Você vai encontrar dois arquivos por cada amostra o r1 e o r2.
Por enquanto preocupe-se com o parâmetro **Per base sequence quality** que dever ser **Q >30**
Você vai ver esse gráfico de barras 📊. Lembra no tutorial de tipos de arquivos que tem a explicação do que é *Phred score*? para refrescar a memória, 🔗[vai lá](https://) **é importante**:exclamation:.
O r2 sempre tem uma qualidade mais baixa.
Aqui neste reporte você também pode confirmar o número de sequências (**Basic Statistics**) que foram analisadas, que deve corresponder à saída do cutadapt
Se quiser conhecer mais sobre este reporte, saber que significam os demais parâmetros. 🔗[Vai aqui](https://dnacore.missouri.edu/PDF/FastQC_Manual.pdf).
### 3 Filtragem das *reads* com baixa qualidade
```coffeescript=
## Crie um novo diretório para armazenar as sequências limpas
mkdir 03.CleanData
## Crie outra pasta onde vai armazenar o que NÃO precisa, o lixo
mkdir unpaired
```
Para a trimagem ✂️ das *reads* com baixa qualidade vai usar o programa Trimmomatic
A continuação alguns parâmetros:
**LEADING** Remove as bases com baixa qualidade do começo da sequência. Especifique a qualidade minima para conservar a base.
**TRAILING** Remove as bases com baixa qualidade do final da sequência. Especifique a qualidade minima para conservar a base.
**SLIDINGWINDOW** Cria uma janela deslizante de determinado número de bases. Elimina a janela completa se a média de qualidade é menor da estipulada. WindowSize:requiredQuality.
**MAXINFO** Executa um aro de qualidade adaptável, equilibrando os benefícios de reter as *reads* com o custo de reter bases com erros. targetLength:strictness.
**MINLEN** Remove *reads* com comprimento menor a o estipulado.
Para mais detalhes consulte o 🔗[manual do Trimmomatic](http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf).
```coffeescript=
## Com este comando você vai processar todas as amostras
for i in 01.PrimersTrim/*1.fq.gz
do
BASE=$(basename $i 1.fq.gz)
## The number of threads is depending of your pc
trimmomatic PE -threads 12 $i 01.PrimersTrim/${BASE}2.fq.gz 03.CleanData/${BASE}1_paired.fq.gz unpaired/${BASE}1_unpaired.fq.gz 03.CleanData/${BASE}2_paired.fq.gz unpaired/${BASE}2_unpaired.fq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MAXINFO:80:0.5 MINLEN:150
done
```
**Atenção!**:exclamation: O parêmetros deste comando podem ser totalmente modificados em função de seus dados e do tipo de sequenciamento (amplicon ou *shotgun*)
Você pode reparar que nas quatro amostras sobreviveram aproximadamente 80% das sequências. Você tem que colocar na balança, poucas reads com boa qualidade ou muitas reads com baixa qualidade. Uma boa porcentagem de sequências para ser levado pra frente é no minimo 60% 👍🏼, melhor se fosse mais de 70% 👌🏼.
Mas... vamos a ver 👀 de novo a qualidade como ficou com o fastqc
```coffeescript=
fastqc -t 10 03.CleanData/* -o 02.FastqcReports/
multiqc 02.FastqcReports/*paired* -o 02.FastqcReports/
## Fazer download desde Filezilla
```
Melhorou né? :thumbsup:
O r2 também melhorou, no entanto não todas estão na caixa verde, mas é aceptável assim
Olhe na seção **Basic Statitics** o número de sequências restantes para atualizar sua planilha.
Até aqui aprendimos como ver 👀 gráficamente a qualidade das *reads* e a filtrar aquelas com baixa qualidade, para levar as melhores para os análises *downstream*.
# PARTE III Sequenciamento *Shotgun*
## Vamos lá! :beginner:
### Ingressar ao servidor. Instruções 🔗[aqui](https://)
Estando na sua 📁 `/data/usuário/` vamos criar uma nova 📁 que vai chamar `shotgun_tutorial`
**Entrada**
```coffeescript=
## Confira onde você está
pwd
```
**Saída**
```coffeescript=
/data/usuário
```
**Entrada**
```coffeescript=
## Crie uma nova pasta
mkdir shotgun_tutorial
## confira
ls
```
**Saída**
```coffeescript=
shotgun_tutorial
```
**Entrada**
```coffeescript=
## Entre na nova pasta
cd shotgun_tutorial/
```
**Todos os comandos vão ser rodados estando dentro desta 📁**
Agora vamos a criar uma 📁 para colocar os dados brutos
```coffeescript=
mkdir 00.RawData
```
### Example dataset
O *example dataset* se encontra numa 📁 no servidor. Vamos a criar um 🔗**link** do dataset exemplo na 📁 que você criou. Imagine que é como um atalho dos que se criam no Windows. Assim, evita-se informação repetida e uso de espaço desnecessariamente.
**Entrada**
```coffeescript=
## Cópiar
ln -s /home/bioinfo/Documentos/examples_datasets/shotgun/* 00.RawData
## Lembre-se que o simbolo * siginifica qualquer caracter. Nesta linha vamos cópiar todos os arquivos dentro da pasta ~/examples_datasets/shotgun/ para a pasta ~/00.RawData/
## Confira
ls 00.RawData/
```
**Saída**
```coffeescript=
amostra1_1.fq amostra1_2.fq amostra2_1.fq amostra2_2.fq amostra3_1.fq amostra3_2.fq
```
**Nota:** Você poderia trabalhar diretamente com seus dados se você tiver, mas a recomendação 🙂 é primeiro praticar com este *dataset*.
### 1 Explorar a qualidade das sequências
```coffeescript=
## Crie um novo diretório para armazenar os reportes de qualidade
mkdir 01.FastqcReports
```
Para analisar 🤓 a qualidade das sequências vamos usar um programa chamado Fastqc, ele faz um scan de cada sequência e determina a qualidade delas segundo vários parâmetros, criando um reporte com interfaz gráfica.
```coffeescript=
## Fastqc está instalado no ambiente virtual bioinfo
source /opt/Miniconda3/bin/activate bioinfo
## Rodar fastqc
fastqc -t 10 00.RawData/* -o 01.FastqcReports/
```
o flag 🚩 -t significa quantos núcleos ou threads você quer que o programa use. Nosso servidor 🖥 tem 22. 10 é suficiente para este processo.
Na anterior linha de comando você processou todas as amostras de uma vez só, simbolizado no caracter *. Lea o comando assim: fastqc por favor me mostra a qualidade das sequências de todas as amostras que se encontram dentro da 📁 `00.RawData/` e coloque os reportes na 📁 `01.FastqcReports/`
Opcionalmente, você pode gerar um reporte único com todas as amostras, usando a ferramenta multiqc
```coffeescript=
multiqc 01.FastqcReports/* -o 01.FastqcReports/
```
Agora para conseguir observar os reportes vai no filezilla e faça download 🔽 deles no seu 💻. Só descarregue os arquivos `.html`
Você vai encontrar dois arquivos por cada amostra o r1 e o r2.
Por enquanto preocupe-se com o parâmetro **Per base sequence quality** que dever ser **Q >30**
Você vai ver esse gráfico de barras 📊. Lembra no tutorial de tipos de arquivos que tem a explicação do que é *Phred score*? para refrescar a memória, 🔗[vai lá](https://) **é importante**:exclamation:.
O r2 sempre tem uma qualidade um pouco mais baixa.
No caso destas amostras, a qualidade delas está muito boa, por tanto não é necessário trimar elas.
Aqui neste reporte você também pode confirmar o número de sequências (**Basic Statistics**) que foram analisadas. Para estas amostras exemplo foram criados subsample de um sequenciamento *shotgun* com 10.000 *reads* para cada amostra. Para facilitar e agilizar as análises.
Se quiser conhecer mais sobre este reporte, saber que significam os demais parâmetros. 🔗[Vai aqui](https://dnacore.missouri.edu/PDF/FastQC_Manual.pdf).
### 2 Filtragem das *reads* com baixa qualidade
Embora as amostras exemplo para esta parte do tutorial tem uma qualidade aceptável para continuar com as análises, vamos fazer o tutorial completo para você aprender os comando para filtragem de sequencias com baixa qualidade.
```coffeescript=
## Crie um novo diretório para armazenar as sequências limpas
mkdir 02.CleanData
## Crie outra pasta onde vai armazenar o que NÃO precisa, o lixo
mkdir unpaired
```
Para a trimagem ✂️ das *reads* com baixa qualidade vai usar o programa Trimmomatic
A continuação alguns parâmetros:
**LEADING** Remove as bases com baixa qualidade do começo da sequência. Especifique a qualidade minima para conservar a base.
**TRAILING** Remove as bases com baixa qualidade do final da sequência. Especifique a qualidade minima para conservar a base.
**SLIDINGWINDOW** Cria uma janela deslizante de determinado número de bases. Elimina a janela completa se a média de qualidade é menor da estipulada. WindowSize:requiredQuality.
**MAXINFO** Executa um aro de qualidade adaptável, equilibrando os benefícios de reter as *reads* com o custo de reter bases com erros. targetLength:strictness.
**MINLEN** Remove *reads* com comprimento menor a o estipulado.
Para mais detalhes consulte o 🔗[manual do Trimmomatic](http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf).
```coffeescript=
## Com este comando você vai processar todas as amostras
for i in 00.RawData/*1.fq
do
BASE=$(basename $i 1.fq)
trimmomatic PE -threads 10 $i 00.RawData/${BASE}2.fq 02.CleanData/${BASE}1_paired.fq unpaired/${BASE}1_unpaired.fq 02.CleanData/${BASE}2_paired.fq unpaired/${BASE}2_unpaired.fq LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MAXINFO:80:0.5 MINLEN:50
done
```
**Atenção!**:exclamation: O parêmetros deste comando podem ser totalmente modificados em função de seus dados e do tipo de sequenciamento (amplicon ou *shotgun*)
O comando deu certo 👍🏼 quando aparece assim:
Você pode reparar que nas três amostras sobreviveram aproximadamente **95%** das sequências. Porque como explicado anteriormente as sequencias *raw* já estavam com boa qualidade. De tudo jeito no caso que tenha que trimar muito, você tem que colocar na balança, poucas reads com boa qualidade ou muitas reads com baixa qualidade. Uma boa porcentagem de sequências para ser levado pra frente é minimo de 60% 👍🏼, melhor se fosse mais de 70% 👌🏼.
Mas... vamos a ver 👀 de novo a qualidade como ficou com o fastqc
```coffeescript=
fastqc -t 10 02.CleanData/* -o 01.FastqcReports/
## multiqc
multiqc 01.FastqcReports/*paired* -o 01.FastqcReports/
## Fazer download desde Filezilla
```
Melhorou né? :thumbsup:
O r2 agora ficou perfeito igual ao r1
Olhe na seção **Basic Statitics** o número de sequências restantes para atualizar sua planilha.
Até aqui aprendimos como ver 👀 gráficamente a qualidade das *reads* e a filtrar aquelas com baixa qualidade, para levar as melhores para os análises *downstream*.
## Referências
📄 **Trimmomatic:** Bolger, A.M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30.
📄 **FastQC:** Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data.
📄 **MultiQC:** Ewels, P., Magnusson, M., Lundin, S., and Käller, M. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047-3048.
**FIM** :sparkle:
Sign in with Wallet
Connect another wallet