# Tutorial para Montagem, Avaliação e Anotação de Genomas
## Tabela de conteúdo
[TOC]
---
A seguir estão detalhadas as intruções para a prática de montagem de genomas. As instruções a seguir apresentam uma sequência de passos para a criação de um projeto de montagem de um genoma. Esta mesma sequência de passos poderá ser adaptada para outros projetos.
## Requisitos
Para um melhor aproveitamento, os seguintes programas devem estar instalados nas estações de trabalho:
* [WinSCP]() - *SCP for Windows*
* [Integrative Genomics Viewer](https://software.broadinstitute.org/software/igv/) - *Genome browser*
* [BRIG](http://brig.sourceforge.net/) - *BLAST Ring Image Generator*
* [BLAST](ftp://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/) - *BLAST*
* *Preferences* > *BRIG Options* > *BLAST binary folder* > *Save settings as default*
* *Preferences* > *BRIG Options* > *Image settings tab* > *Image render memory*
* [Artemis](https://www.sanger.ac.uk/science/tools/artemis) - *Genome browser and annotation tool*
* [AliView](https://ormbunkar.se/aliview/) - *Alignment viewer and editor*
---
## Slides
{%slideshare dgpinheiro/montagem-de-genomas %}
{%slideshare dgpinheiro/predio-gnica %}
{%slideshare dgpinheiro/anotao-gnica-funcional %}
---
## Montagem de genoma Procarioto
Antes da execução da prática de montagem realizaremos as seguintes tarefas:
* [Preparação](#Preparação)
Para os que desejam realizar uma nova simulação:
* [Simulação](#Simulação)
Para os impacientes, seguir para a seguinte etapa:
* [Impacientes](#Impacientes) - Obtendo os dados simulados no formato .fastq
### Preparação
1. Selecionar o diretório $HOME:
```bash=
cd ~
```
2. Criar diretório de análise: montagem do genoma de Candidatus *Nasuia deltocephalinicola* str. NAS-ALF (tamanho do genoma: 112,091 kbps) e entrar no diretório:
```bash=
mkdir ./Ndel
cd ./Ndel
```
3. Criar subdiretórios:
```bash=
mkdir ./raw
mkdir ./ref
```
### Simulação
1. Entrar em ./ref:
```bash=
cd ./ref
```
2. Obter cópia da sequência genômica completa (.fa) e do conjunto de proteínas (.faa):
```bash=
esearch -db nucleotide -query NC_021919.1 | efetch \
-format fasta > Ndel_NC_021919.1.fa
esearch -db nucleotide -query NC_021919.1 | efetch \
-format fasta_cds_aa > Ndel_NC_021919.1.faa
```
3. Criar arquivo de abundância (100% para NC_021919.1)
```bash=
echo -e "NC_021919.1\t1" > abundance.txt
```
4. Gerar "25000" sequências de fragmentos aleatórios em média de "300 pb" com desvio padrão de "30 pb" a partir da sequência genômica "Ndel_NC_021919.1.fa", considerando a abundância contida em "abundance.txt"
```bash=
generate_fragments.py -r ./Ndel_NC_021919.1.fa \
-a ./abundance.txt \
-o ./tmp.frags \
-t 25000 \
-i 300 \
-s 30
```
5. Renomear as sequências para "SAMPLEA" juntamente com um número hexadecimal sequencial começando por "0000000001" e posicionar todas as bases em uma única linha (até 1000 pb). A linha de SED irá substituir a descrição adicional.
```bash=
cat ./tmp.frags.1.fasta | renameSeqs.pl \
-if FASTA \
-of FASTA \
-p SAMPLEA \
-w 1000 | \
sed 's/^>\(\S\+\).*/>\1/' \
> ./frags_A.fa
```
6. Fazer a simulação considerando leituras "paired" 2x "151" bp , com parâmetros configurados ([https://www.ebi.ac.uk/goldman-srv/simNGS/runfiles5/151cycPE/s_4_0099.runfile]) usando simNGS:
```bash=
cat ./frags_A.fa | simNGS -a \
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
-p paired \
/usr/local/bioinfo/simNGS/data/s_4_0099.runfile \
-n 151 > ./SAMPLEA.fastq 2> SAMPLEA.err.txt
```
7. Desintercalar as leituras em dois arquivos SAMPLEA_R1.fastq e SAMPLEA_R2.fastq no diretório ../raw:
```bash=
deinterleave_pairs SAMPLEA.fastq \
-o ../raw/SAMPLEA_R1.fastq \
../raw/SAMPLEA_R2.fastq
```
8. Remover arquivos desnecessários:
```bash=
rm -f ./tmp.frags.1.fasta ./frags_A.fa ./SAMPLEA.fastq ./SAMPLEA.err.txt
```
9. Voltar ao ponto de partida (diretório anterior):
```bash=
cd ../
```
### Impacientes
1. Copiar os arquivos "SAMPLEA_R[12].fastq" a partir do diretório /data/:
```bash=
cp /data/SAMPLEA_R[12].fastq ./raw/
```
## Montagem
### Newbler
Para montagem com Newbler, somente considerar leituras PE no formato /1 e /2 (estilo antigo). Ver explicação [aqui](https://contig.wordpress.com/2011/09/01/newbler-input-iii-a-quick-fix-for-the-new-illumina-fastq-header/).
```bash=
newAssembly -force newbler_asm
addRun -lib PE newbler_asm ./raw/SAMPLEA_R[12].fastq
runProject -mi 90 -ml 30 -a 300 -cpu 1 -scaffold newbler_asm \
> ./newbler_asm/newbler_asm.log.out.txt \
2> ./newbler_asm/newbler_asm.log.err.txt
```
Caso deseje alterar os parâmetros pré-definidos pode consultar este [Manual](http://gensoft.pasteur.fr/docs/454_DataAnalysis/2.9/USM-00058.09_454SeqSys_SWManual-v2.9_PartC.pdf), pg. 27-31.
### Velvet
Utilizaremos no caso do Velvet o script VelvetOptimiser.pl para otimizar dois parâmetros da execução, o cov_cutoff e o tamanho de k. A otimização será a pré-configurada, ou seja,
para a escolha do k-mer 'n50' (métrica N50) - de 27 a 57 com um passo de 16, e para o cov_cutoff 'Lbp' (Total de pares de base em contigs maiores que 1 kbp).
Os parâmetros extras do velvetg são para realizar o passo de scaffolding e utilizar a estimativa de expectativa de cobertura (*).
(*) Ajuda o Velvet a determinar se a leitura é de uma região repetitiva ou única do genoma, ou seja, se a cobertura real é maior do que o esperado e Velvet encontra variações, deverá montar outras versões, caso contrário, se a cobertura é menor do que o esperado, Velvet pode considerar que são leituras do mesmo local do genoma (erros de sequenciamento ou variações polimórficas). Se é menor do que o esperado, e houver uma diferença de base. Ver este [link](https://homolog.us/blogs/genome/2012/06/08/an-explanation-of-velvet-parameter-exp_cov/) para mais explicações.
```bash=
VelvetOptimiser.pl --threads=1 \
--hashs=27 --hashe=57 --step=16 \
--velvetgoptions='-scaffolding yes -exp_cov auto' \
--velvethfiles='-fastq -longPaired raw/SAMPLEA_R1.fastq raw/SAMPLEA_R2.fastq' \
-dir_final ./velvet_asm
```
Caso não consiga uma otimização para *-cov_cutoff* faça a execução do velvetg diretamente no diretório *./velvet_asm* sem utilizar esse parâmetro:
```bash=
velvetg ./velvet_asm -scaffolding yes -exp_cov auto
```
Posteriormente, podemos limpar os demais diretórios e arquivos de Log para as montagens descartadas:
```bash=
rm -fr *_Logfile.txt ./auto_data_*
```
### SPAdes
Nesta execução do SPAdes, vamos deixar de lado a opção de correção de erros (--only-assembler) devido a problemas com este processo nesta versão do SPAdes (v.3.13.1) com este conjunto de dados. Consultar [aqui](https://github.com/ablab/spades/issues/152) para obter mais informações e acompanhar o desenvolvimento do SPAdes. O mesmo process ode montagem com 50 threads de processamento, em outro servidor, demorou 202 minutos.
```bash=
mkdir ./spades_asm_1
spades.py --phred-offset 33 -t 1 --cov-cutoff auto -k auto --only-assembler --careful \
-o ./spades_asm_1 \
-1 ./raw/SAMPLEA_R1.fastq \
-2 ./raw/SAMPLEA_R2.fastq \
> ./spades_asm_1/spades_asm_1.log.out.txt \
2> ./spades_asm_1/spades_asm_1.log.err.txt
```
### MeGAMerge
Montagens fusionadas:
```bash=
MeGAMerge-1.1.pl -cpu=1 \
-single_genome=1 \
-force \
-o=megamerge_asm.fa \
./megamerge_asm \
./spades_asm_1/scaffolds.fasta \
./velvet_asm/contigs.fa \
./newbler_asm/assembly/454Scaffolds.fna
```
## Cobertura
Em profundidade (*Depth coverage*) e em largura (*Breadth coverage*).
<small>Fonte: https://www.danielecook.com/calculate-depth-and-breadth-of-coverage-from-a-bam-file/</small>
## Avaliação
### QUAST
Avaliação das 4 montagens realizadas com QUAST:
1. Criação de links simbólicos
```bash=
ln -s newbler_asm/assembly/454Scaffolds.fna newbler_asm.fa
ln -s velvet_asm/contigs.fa velvet_asm.fa
ln -s spades_asm_1/scaffolds.fasta spades_asm_1.fa
ln -s megamerge_asm/megamerge_asm.fa megamerge_asm.fa
```
2. Execução do QUAST
```bash=
quast.py -r ./ref/Ndel_NC_021919.1.fa --gene-finding \
-1 ./raw/SAMPLEA_R1.fastq \
-2 ./raw/SAMPLEA_R2.fastq \
-m 0 -t 1 \
-o quast_eval \
./newbler_asm.fa \
./velvet_asm.fa \
./spades_asm_1.fa \
./megamerge_asm.fa
```
3. Cópia via SSH dos resultados
```bash=
scp -P 12233 -r <username>@200.17.60.227:~/Ndel/quast_eval/ /tmp/
```
4. Visualização com navegador e acesso ao relatório
Na URL entrar com o caminho:
file:///tmp/quast_eval/report.html
### COMPASS
1. Criar diretório
```bash=
mkdir compass/
```
2. Avaliação com COMPASS (metade da montagem: 56 kbp):
```bash=
compass.pl -n 56000 ./ref/Ndel_NC_021919.1.fa spades_asm_1.fa > compass/spades_asm_1.txt
compass.pl -n 56000 ./ref/Ndel_NC_021919.1.fa velvet_asm.fa > compass/velvet_asm.txt
compass.pl -n 56000 ./ref/Ndel_NC_021919.1.fa newbler_asm.fa > compass/newbler_asm.txt
compass.pl -n 56000 ./ref/Ndel_NC_021919.1.fa megamerge_asm.fa > compass/megamerge_asm.txt
```
3. Visualização de todos os arquivos juntos, somente as 5 linhas (a partir da linha 31) com os valores das métricas:
```bash=
paste ./compass/newbler_asm.txt \
./compass/velvet_asm.txt \
./compass/spades_asm_1.txt \
./compass/megamerge_asm.txt \
| head -n+31 | tail -5 | column -t -s$'\t' | less -S
```
4. Remoção de arquivos desnecessários:
```bash=
rm -f compassRun.*
```
## Predição gênica
A predição gênica será realizada com PROKKA para a montagem selecionada (megamerge_asm.fa)
1. Execução do PROKKA
```bash=
prokka --cpus 1 \
--locustag Ndel \
--rnammer \
--kingdom Bacteria \
--gram neg \
--outdir prokka_pred \
--force --prefix \
Ndel-1.0 megamerge_asm.fa
```
### BUSCO
Para avaliação da completude, utilizaremos o BUSCO com as proteínas preditas (./prokka_pred/Ndel-1.0.faa).
1. Execução do BUSCO com genes essenciais da linhagem Bacteria
* BUSCO para a referência
```bash=
run_BUSCO.py --cpu 1 --in ./ref/Ndel_NC_021919.1.faa \
--out ref \
-l /usr/local/bioinf/busco_lineages/bacteria_odb9/ \
-m prot --force
```
* BUSCO para a montagem
```bash=
run_BUSCO.py --cpu 1 --in ./prokka_pred/Ndel-1.0.faa \
--out Ndel-1.0 --tmp_path /tmp/ \
-l /usr/local/bioinf/busco_lineages/bacteria_odb9/ \
-m prot --force
```
2. Gerando Gráfico
* Criar diretório
```bash=
mkdir ./busco_eval
```
* Criar links simbólicos
```bash=
ln -s ../run_ref/short_summary_ref.txt busco_eval/short_summary_ref.txt
ln -s ../run_Ndel-1.0/short_summary_Ndel-1.0.txt busco_eval/short_summary_Ndel-1.0.txt
```
* Executar
```bash=
generate_plot.py -wd ./busco_eval/
```
3. Copiar via SSH dos resultados
```bash=
scp -P 12233 -r <username>@200.17.60.227:~/Ndel/busco_eval/ /tmp/
```
4. Visualização com navegador e acesso à imagem
Na URL entrar com o caminho:
file:///tmp/busco_eval/busco_figure.png
5. Avaliação com BLAST
* Ir para o diretório ref/
```bash=
cd ./ref
```
* Criar índice para o conjunto de proteínas preditas no NCBI (referência)
```bash=
makeblastdb -dbtype prot -in Ndel_NC_021919.1.faa
```
* Executar BLASTP predição do prokka X predição do NCBI
* Preparação
```=bash
cd ../
mkdir blast_cmp/
cd blast_cmp/
```
* Execução
```bash=
blastp -db ../ref/Ndel_NC_021919.1.faa \
-query ../prokka_pred/Ndel-1.0.faa \
-num_threads 1 -max_target_seqs 1 -outfmt 6 > blast.txt
```
* Comparação
```bash=
grep '^>' ../ref/Ndel_NC_021919.1.faa | \
sed 's/^>//' | sed 's/ .*//' | \
sort -u > ALL_NCBI.txt
grep '^>' ../prokka_pred/Ndel-1.0.faa | \
sed 's/^>//' | sed 's/ .*//' | \
sort -u > ALL_PROKKA.txt
cut -f 2 blast.txt | sort -u > MATCHES_NCBI.txt
cut -f 1 blast.txt | sort -u > MATCHES_PROKKA.txt
```
Identificadas pelo NCBI e não pelo Prokka
```bash=
diff ALL_NCBI.txt MATCHES_NCBI.txt | grep '^<'
```
Identificadas pelo Prokka e não pelo NCBI
```bash=
diff ALL_PROKKA.txt MATCHES_PROKKA.txt | grep '^<'
```
## Anotação Gênica Funcional
Antes de executar os programas de anotação, podemos separar as sequências de proteínas preditas, em arquivos com até 100 sequências (200 linhas - por isso cada sequência deve conter uma linha de cabeçalho e uma de sequência [-l 200]). Os aquivos resultantes terão um prefixo com autocontagem com 3 caracteres [-a 3]:
1. Criar diretório com as predições em vários arquivos
```bash=
mkdir chunks_annot
```
2. Executar os comandos abaixo para gerar os arquivos com até 100 sequências
```bash=
oneLineFasta.pl -i prokka_pred/Ndel-1.0.faa | split -l 200 -a 3 -d - chunks_annot/input_file.chunk_
```
3. Verificar a contagem de sequências por partição
```bash=
grep -c '^>' chunks_annot/input_file.chunk_*
```
3. InterProScan
* Criar diretório
```bash=
mkdir ./interproscan_annot/
```
* Executar
```bash=
for f in ./chunks_annot/*.chunk_*; do
bn=`basename ${f}`;
interproscan.sh --goterms --iprlookup --pathways --seqtype p -f TSV \
--input ${f} \
--outfile interproscan_annot/${bn}.iprscan.tsv;
done
```
* Concatenar os resultados para Ndel-1.0-iprscan.tsv
```bash=
cat interproscan_annot/*.iprscan.tsv > ./interproscan_annot/Ndel-1.0-iprscan.tsv
```
* Retirar do resultado somente a informação de ID, IPR e descrição para um outro arquivo
```bash=
cut -f 1,12,13 interproscan_annot/Ndel-1.0-iprscan.tsv | \
grep -P '\tIPR[0-9]+\t' |\
nsort -u > interproscan_annot/Ndel-1.0-iprscan-IPR.tsv
```
* Agregar IPRs e descrições de uma mesma proteína em uma única linha
```bash=
aggregatoR.R --x=interproscan_annot/Ndel-1.0-iprscan-IPR.tsv \
--sep=";" \
--colnames.x='id,ipr,ipr.desc' \
--by.x='id' \
--noh.x \
--out=interproscan_annot/Ndel-1.0-iprscan-IPR-agg.tsv
```
2. eggNOG-mapper
* Criar diretório
```bash=
mkdir ./emapper_annot
```
* Executar
* Phase 1. Homology searches*
```bash=
for f in ./chunks_annot/*.chunk_*; do
bn=`basename ${f}`;
emapper.py -m diamond --no_annot --no_file_comments --cpu 1 -i ${f} -o ./emapper_annot/${bn};
done
```
* Concatenar os resultados para Ndel-1.0-iprscan.tsv
```bash=
cat emapper_annot/*.seed_orthologs > emapper_annot/Ndel-1.0-emapper.seed_orthologs
```
* Executar
* Phase 2. Orthology and functional* annotation
```bash=
emapper.py --cpu 1 --annotate_hits_table emapper_annot/Ndel-1.0-emapper.seed_orthologs \
--no_file_comments \
-o emapper_annot/Ndel-1.0-emapper.tsv
```
* Retirar do resultado somente a informação de ID, símbolo do gene e descrição para um outro arquivo
```bash=
cut -f 1,6,22 emapper_annot/Ndel-1.0-emapper.tsv.emapper.annotations > emapper_annot/Ndel-1.0-emapper.tsv.emapper.annotations-desc.txt
```
* Agregar símbolos gênicos e descrições de uma mesma proteína em uma única linha
```bash=
aggregatoR.R --x=emapper_annot/Ndel-1.0-emapper.tsv.emapper.annotations-desc.txt \
--sep=";" \
--colnames.x='id,eggnog.gene,eggnog.desc' \
--by.x='id' \
--noh.x \
--out=emapper_annot/Ndel-1.0-emapper.tsv.emapper.annotations-desc-agg.txt
```
3. Reunir as principais descrições em um único arquivo
* Criar diretório
```bash=
mkdir ./merged_annot
```
* Primeira fusão: PROKKA + InterProScan
```bash=
mergeR.R --x=./prokka_pred/Ndel-1.0.tsv \
--y=./interproscan_annot/Ndel-1.0-iprscan-IPR-agg.tsv \
--by.x='locus_tag' \
--by.y='id' \
--all.x \
--out=./merged_annot/prokka_interproscan.txt
```
* Segunda fusão: (PROKKA + InterProScan) + eggNOG-mapper
```bash=
mergeR.R --x=./merged_annot/prokka_interproscan.txt \
--y=./emapper_annot/Ndel-1.0-emapper.tsv.emapper.annotations-desc-agg.txt \
--by.x='locus_tag' \
--by.y='id' \
--all.x \
--out=./merged_annot/prokka_interproscan_eggnog.txt
```
* Consultar o resultado final da anotação (visual tabular):
```bash=
cat ./merged_annot/prokka_interproscan_eggnog.txt | column -t -s$'\t' | less -S
```
* Transformar para .xls (Excel)
```bash=
tsv2xls.pl -i ./merged_annot/prokka_interproscan_eggnog.txt -o ./merged_annot/prokka_interproscan_eggnog -h
```
* Cópia via SSH dos resultados
```bash=
scp -P 12233 -r <username>@200.17.60.227:~/Ndel/merged_annot/ /tmp/
```
> [name=Daniel Guariz Pinheiro] Encerramos aqui a prática de Montagem, Predição e Anotação realizadas a partir de um conjunto de dados gerados por simulação de sequenciamento de uma bactéria endosimbionte (Candidatus *Nasuia deltocephalinicola* str. NAS-ALF).
---
## Montagem de genoma Eucarioto
### Preparação
1. Selecionar o diretório $HOME:
```bash=
cd ~
```
2. Criar diretório de análise: montagem do cromossomo III do genoma de *Schizosaccharomyces pombe* (*fission yeast*) (tamanho do cromossomo: 2,45Mb) e entrar no diretório:
```bash=
mkdir ./Spom
cd ./Spom
```
3. Criar subdiretórios:
```bash=
mkdir ./raw
mkdir ./ref
```
### Simulação
1. Entrar em ./ref:
```bash=
cd ./ref
```
2. Obter cópia da sequência genômica completa (.fa) e do conjunto de proteínas (.faa):
```bash=
esearch -db nucleotide -query NC_003421.2 | efetch \
-format fasta > Spom_NC_003421.2.fa
esearch -db nucleotide -query NC_003421.2 | efetch \
-format fasta_cds_aa > Spom_NC_003421.2.faa
```
3. Criar arquivo de abundância (100% para NC_003421.2)
```bash=
echo -e "NC_003421.2\t1" > abundance.txt
```
4. Gerar "100000" sequências de fragmentos aleatórios em média de "300 pb" com desvio padrão de "30 pb" a partir da sequência genômica "Spom_NC_003421.2.fa", considerando a abundância contida em "abundance.txt"
```bash=
generate_fragments.py -r ./Spom_NC_003421.2.fa \
-a ./abundance.txt \
-o ./tmp.frags \
-t 100000 \
-i 300 \
-s 30
```
5. Renomear as sequências para "SAMPLEB" juntamente com um número hexadecimal sequencial começando por "0000000001" e posicionar todas as bases em uma única linha (até 1000 pb). A linha de SED irá substituir a descrição adicional.
```bash=
cat ./tmp.frags.1.fasta | renameSeqs.pl \
-if FASTA \
-of FASTA \
-p SAMPLEB \
-w 1000 | \
sed 's/^>\(\S\+\).*/>\1/' \
> ./frags_B.fa
```
6. Fazer a simulação considerando leituras "paired" 2x "151" bp , com parâmetros configurados ([https://www.ebi.ac.uk/goldman-srv/simNGS/runfiles5/151cycPE/s_4_0099.runfile]) usando simNGS:
```bash=
cat ./frags_B.fa | simNGS -a \
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
-p paired \
/usr/local/bioinfo/simNGS/data/s_4_0099.runfile \
-n 151 > ./SAMPLEB.fastq 2> SAMPLEB.err.txt
```
7. Desintercalar as leituras em dois arquivos SAMPLEB_R1.fastq e SAMPLEB_R2.fastq no diretório ../raw:
```bash=
deinterleave_pairs SAMPLEB.fastq \
-o ../raw/SAMPLEB_R1.fastq \
../raw/SAMPLEB_R2.fastq
```
8. Remover arquivos desnecessários:
```bash=
rm -f ./tmp.frags.1.fasta ./frags_B.fa ./SAMPLEB.fastq ./SAMPLEB.err.txt
```
9. Voltar ao ponto de partida (diretório anterior):
```bash=
cd ../
```
### Impacientes
1. Copiar os arquivos "SAMPLEA_R[12].fastq" a partir do diretório /data/:
```bash=
cp /data/SAMPLEB_R[12].fastq ./raw/
```
## Montagem
### Newbler
Para montagem com Newbler, somente considerar leituras PE no formato /1 e /2 (estilo antigo). Ver explicação [aqui](https://contig.wordpress.com/2011/09/01/newbler-input-iii-a-quick-fix-for-the-new-illumina-fastq-header/).
```bash=
newAssembly -force newbler_asm
addRun -lib PE newbler_asm ./raw/SAMPLEB_R[12].fastq
runProject -mi 90 -ml 30 -a 300 -cpu 1 -scaffold newbler_asm \
> ./newbler_asm/newbler_asm.log.out.txt \
2> ./newbler_asm/newbler_asm.log.err.txt
```
Caso deseje alterar os parâmetros pré-definidos pode consultar este [Manual](http://gensoft.pasteur.fr/docs/454_DataAnalysis/2.9/USM-00058.09_454SeqSys_SWManual-v2.9_PartC.pdf), pg. 27-31.
### Velvet
Utilizaremos no caso do Velvet o script VelvetOptimiser.pl para otimizar dois parâmetros da execução, o cov_cutoff e o tamanho de k. A otimização será a pré-configurada, ou seja,
para a escolha do k-mer 'n50' (métrica N50) - de 27 a 57 com um passo de 16, e para o cov_cutoff 'Lbp' (Total de pares de base em contigs maiores que 1 kbp).
Os parâmetros extras do velvetg são para realizar o passo de scaffolding e utilizar a estimativa de expectativa de cobertura (*).
(*) Ajuda o Velvet a determinar se a leitura é de uma região repetitiva ou única do genoma, ou seja, se a cobertura real é maior do que o esperado e Velvet encontra variações, deverá montar outras versões, caso contrário, se a cobertura é menor do que o esperado, Velvet pode considerar que são leituras do mesmo local do genoma (erros de sequenciamento ou variações polimórficas). Se é menor do que o esperado, e houver uma diferença de base. Ver este [link](https://homolog.us/blogs/genome/2012/06/08/an-explanation-of-velvet-parameter-exp_cov/) para mais explicações.
```bash=
VelvetOptimiser.pl --threads=1 \
--hashs=27 --hashe=57 --step=16 \
--velvetgoptions='-scaffolding yes -exp_cov auto' \
--velvethfiles='-fastq -longPaired raw/SAMPLEB_R1.fastq raw/SAMPLEB_R2.fastq' \
-dir_final ./velvet_asm
```
Caso não consiga uma otimização para *-cov_cutoff* faça a execução do velvetg diretamente no diretório *./velvet_asm* sem utilizar esse parâmetro:
```bash=
velvetg ./velvet_asm -scaffolding yes -exp_cov auto
```
Posteriormente, podemos limpar os demais diretórios e arquivos de Log para as montagens descartadas:
```bash=
rm -fr *_Logfile.txt ./auto_data_*
```
### SPAdes
Nesta execução do SPAdes, vamos deixar de lado a opção de correção de erros (--only-assembler) devido a problemas com este processo nesta versão do SPAdes (v.3.13.1) com este conjunto de dados. Consultar [aqui](https://github.com/ablab/spades/issues/152) para obter mais informações e acompanhar o desenvolvimento do SPAdes. O mesmo process ode montagem com 50 threads de processamento, em outro servidor, demorou 202 minutos.
```bash=
mkdir ./spades_asm_1
spades.py --phred-offset 33 -t 1 --cov-cutoff auto -k auto --only-assembler --careful \
-o ./spades_asm_1 \
-1 ./raw/SAMPLEB_R1.fastq \
-2 ./raw/SAMPLEB_R2.fastq \
> ./spades_asm_1/spades_asm_1.log.out.txt \
> ./spades_asm_1/spades_asm_1.log.err.txt
```
### MeGAMerge
Montagens fusionadas:
```bash=
MeGAMerge-1.1.pl -cpu=1 \
-single_genome=1 \
-force \
-o=megamerge_asm.fa \
./megamerge_asm \
./spades_asm_1/scaffolds.fasta \
./velvet_asm/contigs.fa \
./newbler_asm/assembly/454Scaffolds.fna
```
## Avaliação
### QUAST
Avaliação das 3 montagens realizadas com QUAST:
1. Criação de links simbólicos
```bash=
ln -s newbler_asm/assembly/454Scaffolds.fna newbler_asm.fa
ln -s velvet_asm/contigs.fa velvet_asm.fa
ln -s spades_asm_1/scaffolds.fasta spades_asm_1.fa
ln -s megamerge_asm/megamerge_asm.fa megamerge_asm.fa
```
2. Execução do QUAST
Note os parâmetros --eukaryote e --fungus que influenciarão na predição gênica realizada internamente pelo QUAST com o programa GeneMark-ES(*).
```bash=
quast.py -r ./ref/Spom_NC_003421.2.fa \
--gene-finding \
--eukaryote \
--fungus \
-1 ./raw/SAMPLEB_R1.fastq \
-2 ./raw/SAMPLEB_R2.fastq \
-m 0 -t 1 \
-o quast_eval \
newbler_asm.fa \
velvet_asm.fa \
spades_asm_1.fa \
megamerge_asm.fa
```
(*) A versão original do programa QUAST foi alterada para indicar ao GeneMark-ES para realizar a etapa de treinamento com contigs de no mínimo 500 pb. Caso contrário, ele aponta erros na predição com Velvet.
3. Cópia via SSH dos resultados
```bash=
scp -P 12233 -r <username>@200.17.60.227:~/Ndel/quast_eval/ /tmp/
```
4. Visualização com navegador e acesso ao relatório
Na URL entrar com o caminho:
file:///tmp/quast_eval/report.html
### COMPASS
1. Criar diretório
```bash=
mkdir compass/
```
2. Avaliação com COMPASS (10% da montagem do cromossomo: 255000 bps):
```bash=
compass.pl -n 255000 ./ref/Spom_NC_003421.2.fa spades_asm_1.fa > compass/spades_asm_1.txt
compass.pl -n 255000 ./ref/Spom_NC_003421.2.fa velvet_asm.fa > compass/velvet_asm.txt
compass.pl -n 255000 ./ref/Spom_NC_003421.2.fa newbler_asm.fa > compass/newbler_asm.txt
compass.pl -n 255000 ./ref/Spom_NC_003421.2.fa megamerge_asm.fa > compass/megamerge_asm.txt
```
3. Visualização de todos os arquivos juntos, somente as 5 linhas (a partir da linha 31) com os valores das métricas:
```bash=
paste ./compass/newbler_asm.txt \
./compass/velvet_asm.txt \
./compass/spades_asm_1.txt \
./compass/megamerge_asm.txt \
| head -n+31 | tail -5
```
4. Remoção de arquivos desnecessários:
```bash=
rm -f compassRun.*
```
### REAPR
:::info
**[ IMPORTANTE ]**: REAPR aceita apenas *reads* com o mesmo tamanho!
**[ DICA ]:** PRINSEQ *trima* sequências por tamanho.
:::
* Criar diretório
```bash=
mkdir -p ./reapr_eval/prinseq_proc/
mkdir -p ./reapr_eval/bowtie/idx/
mkdir -p ./reapr_eval/FragmentSize/
mkdir -p ./reapr_eval/perfectmap/
mkdir -p ./reapr_eval/smaltmap/
```
* Executar
```bash=
prinseq-lite.pl -fastq ./raw/SAMPLEB_R1.fastq \
-fastq2 ./raw/SAMPLEB_R2.fastq \
-trim_to_len 150 \
-out_format 3 \
-out_bad null \
-out_good ./reapr_eval/prinseq_proc/SAMPLEB \
-min_len 150
```
* Entrar no diretório onde será gerado índice para montagem
```bash=
cd ./reapr_eval/bowtie/idx/
```
* Criar link simbólico do fasta da montagem, neste caso montagem com SPAdes.
```bash=
ln -s ../../../spades_asm_1.fa asm.fa
```
* Indexar a montagem
```bash=
bowtie2-build asm.fa asm
```
* Trocar de diretório para reapr_eval/
```bash=
cd ../../
```
* Mapeamento das *reads* na montagem
```bash=
bowtie2 -p 1 -x ./bowtie/idx/asm \
-1 ./prinseq_proc/SAMPLEB_1.fastq \
-2 ./prinseq_proc/SAMPLEB_2.fastq \
-S ./bowtie/mapped.sam \
&> ./bowtie/bowtie.out
```
* Conversão de formato SAM >> BAM
```bash=
samtools view -@ 1 -S -b bowtie/mapped.sam > bowtie/mapped.bam
```
* Ordenando o arquivo BAM e indexando
```bash=
samtools sort -@ 1 bowtie/mapped.bam > bowtie/mapped.sorted.bam
samtools index bowtie/mapped.sorted.bam
```
* Avaliação do tamanho do fragmento: deeptools - bamPEFragmentSize
```bash=
bamPEFragmentSize -hist FragmentSize/fragmentSize.png \
-T "Fragment size of PE data" \
--maxFragmentLength 500 \
-b bowtie/mapped.sorted.bam \
&> FragmentSize/bamPEFragmentSize.out
```
**Resultado (FragmentSize.png)**
Checar o tamanho médio do fragmento para usar nas etapas anteriores (REAPR).
* REAPR - Ferramenta universal para avaliação de montagens de genoma
- [**ARTIGO:**](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2013-14-5-r47) *REAPR: a universal tool for genome assembly evaluation*
- [**MANUAL:**](ftp://ftp.sanger.ac.uk/pub/resources/software/reapr/Reapr_1.0.17.manual.pdf) REAPR version 1.0.17
-- Link simbólico
```bash=
ln -s ../spades_asm_1.fa asm.fa
```
-- FACHECK
Checa os nomes das sequência (cabeçalho do arquivos fasta) para não interromper a etapa da pipeline.
```bash=
reapr facheck asm.fa
```
:::info
Se retornar um erro, faça um novo arquivo FASTA da sua montagem editando o cabeçalho. Alguns caracteres especiais não são aceitos.
:::
-- PERFECTMAP
Alinhamentos perfeitos com a montagem.
```bash=
reapr perfectmap asm.fa \
~/Spom/reapr_eval/prinseq_proc/SAMPLEB_1.fastq \
~/Spom/reapr_eval/prinseq_proc/SAMPLEB_2.fastq \
300 \
perfectmap/perfect &> perfectmap/perfect.out
```
:::warning
Os caminhos dos arquivos ".fastq" acima devem ser absolutos e não relativos.
:::
-- SMALTMAP
Alinhamentos gerais com a montagem, não somente as melhores.
```bash=
reapr smaltmap -n 1 -u 300 asm.fa \
~/Spom/reapr_eval/prinseq_proc/SAMPLEB_1.fastq \
~/Spom/reapr_eval/prinseq_proc/SAMPLEB_2.fastq \
smaltmap/smalt.mapped.bam &> smaltmap/smaltmap.out
```
-- PIPELINE
Avaliação dos alinhamentos e pontos de quebra (*breakpoints*). Depois da avaliação fazer a quebra ( -break b=1 ).
```bash=
reapr pipeline -break b=1 asm.fa smaltmap/smalt.mapped.bam pipeline perfectmap/perfect
```
:::info
Arquivo com a montagem quebrada: **./pipeline/04.break.broken_assembly.fa**
Arquivo com a montagem original:
**./00.assembly.fa**
Arquivo GFF com os erros:
**./pipeline/03.score.errors.gff.gz**
Arquivo BAM com o alinhamento:
**./pipeline/00.in.bam**
Arquivo BAI com o índice do alinhamento
**./pipeline/00.in.bam.bai**
:::
(*) Os arquivos podem ser carregados no IGV para visualização
### SSPACE
[**ARTIGO:**](https://academic.oup.com/bioinformatics/article/27/4/578/197626) *Scaffolding pre-assembled contigs using SSPACE*
[**GITHUB**](https://github.com/nsoranzo/sspace_basic)
- Entrar no diretório para criação do diretório
```bash=
cd ~/Spom/
mkdir ./sspace_asm
cd ./sspace_asm
```
- Consultar informações do fragmento para criar o arquivo com informações da biblioteca
```bash=
more ../reapr_eval/FragmentSize/bamPEFragmentSize.out
```
:::info
Fragment lengths:
Min.: 0.0
1st Qu.: 268.0
Mean: 266.278061224
Median: 295.0
3rd Qu.: 317.0
Max.: 392.0
Std: 98.1773526834
MAD: 24.0
:::
- Criar um arquivo com contendo informações das bibliotecas
```bash=
echo "lib1 $(echo ~/Spom/raw/SAMPLEB_R1.fastq) $(echo ~/Spom/raw/SAMPLEB_R2.fastq) 266 0.25 FR" > ./libraries.txt
```
:::warning
**Colunas**
--Nome da biblioteca
--Arquivo_R1 (Caminho do arquivo das reads *forward*)
--Arquivo_R2 (Caminho do arquivo das reads *reverse*)
--Tamanho do inserto
--Proporção de erro em relação ao tamanho do fragmento
--Orientação
:::
```bash=
SSPACE_Basic.pl \
-s ../reapr_eval/pipeline/04.break.broken_assembly.fa \
-l libraries.txt \
-x 1 -r 0.6 -k 2 -o 1 -n 10 -z 300
```
:::warning
-x 1 (*extension using *paired-end reads*)
-r 0.6 (*minimum base ratio used to accept a overhang consensus base*)
-k 1 (*minimum number of links - read pairs - to compute scaffold*)
-o 1 (*minimum number of reads needed to call a base during an extension*)
-n 5 (*minimum overlap required between contigs to merge adjacent contigs in a scaffold*)
-z (*minimum contig length used for scaffolding. Filters out contigs below this value*)
:::
:::info
Arquivo com a montagem final: **standard_output.final.scaffolds.fasta**
:::
## Predição gênica
A predição gênica será realizada com BRAKER2 para a montagem selecionada (standard_output.final.scaffolds.fasta)
1. Trocar de diretório e criar diretório
```bash=
cd ../
mkdir ./braker2_pred
cd ./braker2_pred
```
3. Execução do BRAKER2
BRAKER2 é uma extensão do BRAKER1 o qual permite uma automatização do treinamento de preditores gênicos GeneMark-EX e AUGUSTUS a patir de dados de RNA-Seq (não é o nosso caso) e/ou informação de homologia (utilizaremos aqui o conjunto de proteínas da mesma espécie), e capaz de integrar essas evidências na predição.
- Treinamento / Predição com BRAKER2
:::danger
NÃO EXECUTAR POIS DEMORA MUITO! IR PARA PREDIÇÃO COM AUGUSTUS.
:::
```bash=
braker.pl --species=Spom \
--genome=../sspace_asm/standard_output.final.scaffolds.fasta \
--prot_seq=../ref/Spom_NC_003421.2.faa \
--prg=gth --trainFromGth --gth2traingenes \
--fungus \
--softmasking
```
:::danger
Se já treinado (como já foi o caso!), ou se deseja utilizar uma espécie já treinada (como iremos utilizar abaixo - saccharomyces_cerevisiae_S288C) com augustus diretamente.
:::
- Predição com Augustus
```bash=
augustus --gff3=on \
--species=saccharomyces_cerevisiae_S288C \
../sspace_asm/standard_output.final.scaffolds.fasta \
2> ./augustus.err.txt > augustus.gff
```
Depois do treinamento do modelo gênico para Spom com BRAKER2, é possível fazer a predição usando este novo modelo, a partir do augustus
```bash=
augustus --gff3=on \
--species=Spom \
../sspace_asm/standard_output.final.scaffolds.fasta \
2> ./Spom_augustus.err.txt > Spom_augustus.gff
```
:::info
Para encontrar as espécies já treinadas no Augustus:
```bash=
augustus --species=help
```
:::
- Obter FASTA com as proteínas preditas
```bash=
getAnnoFasta.pl augustus.gff
```
### BUSCO
Para avaliação da completude, utilizaremos o BUSCO com as proteínas preditas.
1. Trocar para o diretório base do projeto de montagem
```bash=
cd ~/Spom
```
2. Execução do BUSCO com genes essenciais da linhagem Bacteria
* Para a referência
```bash=
run_BUSCO.py --cpu 1 --in ./ref/Spom_NC_003421.2.faa \
--out ref \
-l /usr/local/bioinf/busco_lineages/ascomycota_odb9 \
-m prot --force
```
* Para a montagem
```bash=
run_BUSCO.py --cpu 1 \
--in ./braker2_pred/augustus.aa \
--out Spom-1.0 \
-l /usr/local/bioinf/busco_lineages/ascomycota_odb9 \
-m prot --force
```
2. Gerando Gráfico
* Criar diretório
```bash=
mkdir ./busco_eval
```
* Criar links simbólicos
```bash=
ln -s ../run_ref/short_summary_ref.txt busco_eval/short_summary_ref.txt
ln -s ../run_Spom-1.0/short_summary_Spom-1.0.txt busco_eval/short_summary_Spom-1.0.txt
```
* Executar
```bash=
generate_plot.py -wd ./busco_eval/
```
3. Copiar via SSH dos resultados
```bash=
scp -P 12233 -r <username>@200.17.60.227:~/Ndel/busco_eval/ /tmp/
```
4. Visualização com navegador e acesso à imagem
Na URL entrar com o caminho:
file:///tmp/busco_eval/busco_figure.png
## Anotação Gênica Funcional
Antes de executar os programas de anotação, podemos separar as sequências de proteínas preditas, em arquivos com até 100 sequências (200 linhas - por isso cada sequência deve conter uma linha de cabeçalho e uma de sequência [-l 200]). Os aquivos resultantes terão um prefixo com autocontagem com 3 caracteres [-a 3]:
1. Criar diretório com as predições em vários arquivos
```bash=
mkdir chunks_annot
```
2. Executar os comandos abaixo para gerar os arquivos com até 100 sequências
```bash=
oneLineFasta.pl -i ./braker2_pred/augustus.aa | split -l 200 -a 3 -d - chunks_annot/input_file.chunk_
```
3. InterProScan
* Criar diretório
```bash=
mkdir ./interproscan_annot/
```
* Executar
```bash=
for f in ./chunks_annot/*.chunk_*; do
bn=`basename ${f}`;
interproscan.sh --goterms --iprlookup --pathways --seqtype p -f TSV \
--input ${f} \
--outfile interproscan_annot/${bn}.iprscan.tsv;
done
```
* Concatenar os resultados para Spom-1.0-iprscan.tsv
```bash=
cat interproscan_annot/*.iprscan.tsv > ./interproscan_annot/Spom-1.0-iprscan.tsv
```
* Retirar do resultado somente a informação de ID, IPR e descrição para um outro arquivo
```bash=
cut -f 1,12,13 interproscan_annot/Spom-1.0-iprscan.tsv | \
grep -P '\tIPR[0-9]+\t' |\
nsort -u > interproscan_annot/Spom-1.0-iprscan-IPR.tsv
```
* Agregar IPRs e descrições de uma mesma proteína em uma única linha
```bash=
aggregatoR.R --x=interproscan_annot/Spom-1.0-iprscan-IPR.tsv \
--sep=";" \
--colnames.x='id,ipr,ipr.desc' \
--by.x='id' \
--noh.x \
--out=interproscan_annot/Spom-1.0-iprscan-IPR-agg.tsv
```
2. eggNOG-mapper
* Criar diretório
```bash=
mkdir ./emapper_annot
```
* Executar
* *Phase 1. Homology searches*
```bash=
for f in ./chunks_annot/*.chunk_*; do
bn=`basename ${f}`;
emapper.py -m diamond --no_annot --no_file_comments --cpu 1 -i ${f} -o ./emapper_annot/${bn};
done
```
* Concatenar os resultados para Ndel-1.0-iprscan.tsv
```bash=
cat emapper_annot/*.seed_orthologs > emapper_annot/Spom-1.0-emapper.seed_orthologs
```
* Executar
* *Phase 2. Orthology and functional* annotation
```bash=
emapper.py --cpu 1 --annotate_hits_table emapper_annot/Spom-1.0-emapper.seed_orthologs \
--no_file_comments \
-o emapper_annot/Spom-1.0-emapper.tsv
```
* Retirar do resultado somente a informação de ID, símbolo do gene e descrição para um outro arquivo
```bash=
cut -f 1,6,22 emapper_annot/Spom-1.0-emapper.tsv.emapper.annotations > emapper_annot/Spom-1.0-emapper.tsv.emapper.annotations-desc.txt
```
* Agregar símbolos gênicos e descrições de uma mesma proteína em uma única linha
```bash=
aggregatoR.R --x=emapper_annot/Spom-1.0-emapper.tsv.emapper.annotations-desc.txt \
--sep=";" \
--colnames.x='id,eggnog.gene,eggnog.desc' \
--by.x='id' \
--noh.x \
--out=emapper_annot/Spom-1.0-emapper.tsv.emapper.annotations-desc-agg.txt
```
3. Reunir as principais descrições em um único arquivo
* Criar diretório
```bash=
mkdir ./merged_annot
```
* Fusão InterProScan + eggNOG-mapper
```bash=
mergeR.R --x=./interproscan_annot/Spom-1.0-iprscan-IPR-agg.tsv \
--y=./emapper_annot/Spom-1.0-emapper.tsv.emapper.annotations-desc-agg.txt \
--by.x='id' \
--by.y='id' \
--all.x \
--all.y \
--out=./merged_annot/interproscan_eggnog.txt
```
* Consultar o resultado final da anotação (visual tabular):
```bash=
cat ./merged_annot/interproscan_eggnog.txt | column -t -s$'\t' | less -S
```
* Transformar para .xls (Excel)
```bash=
tsv2xls.pl -i ./merged_annot/interproscan_eggnog.txt -o ./merged_annot/interproscan_eggnog -h
```
* Cópia via SSH dos resultados
```bash=
scp -P 12233 -r <username>@200.17.60.227:~/Ndel/merged_annot/ /tmp/
```
> [name=Daniel Guariz Pinheiro] Encerramos aqui a prática de Montagem, Predição e Anotação realizadas a partir de um conjunto de dados gerados por simulação de sequenciamento de um cromossomo de levedura (*Schizosaccharomyces pombe*).
Sign in with Wallet
Connect another wallet