Montagem de novo e anotação de genomas microbianos

Olá pessoal! Neste tutorial mostrarei uma pequena pipeline para montagem e anotação de genomas microbianos sequenciados com next generation sequencing (NGS). A ideia é desmistificar o processo de análise de dados de sequências, mostrando como é possível ir de dados brutos a um genoma pronto para ser submetido para o GenBank :D

WGS: Whole Genome Shotgun

Hoje em dia, com o advento do sequenciamento de nova geração (NGS), é possível se analisar um genoma microbiano completo em poucas horas dependendo da plataforma que está sendo utilizada. No caso de genomas pequenos, de organismo procarióticos, a plataforma mais utilizadas é a Illumina MiSeq, que utiliza o método de sequenciamento por síntese (sequencing by synthesis) como base. Esta plataforma permite a geração de leituras com tamanho de ~300 bp e uso de bibliotecas do tipo paired-end (onde cada fragmento de DNA dá origem a 2 leituras, uma para cada extremidade), o que facilita o processo de montagem e permite diversas etapas adicionais para scaffolding e fechamento de gaps.

Após o sequenciamento as milhões de leituras geradas são processadas para se reconstruir a sequência do genoma original. No caso do NGS este processo, denominado montagem (assembly), é geralmente realizado por algoritmos que usam como base os grafos de Bruijn. Esta estrutura de dados é construida a partir de fragmentos de tamanho fixo (k) gerados a partir das leituras de sequenciamento, sendo estes denominados k-mers.

A estrutura de grafo é construída a partir de sobreposições “k-1”. Se o k for igual a 21, por exemplo, todas as sobreposições de 20 base identificadas entre os k-mers são utilizadas na construção de conexões no grafo. Após isso, o mesmo é percorrido de modo a se encontrar os menor conjunto de caminhos independentes.

O processo de montagem resulta na produção de contigs, que são sequencias contiguas, e portanto, ininterruptas. Estas sequências não possuem, por exemplo, regiões desconhecidas (gaps), pois são geradas a partir da sobreposição direta de leituras ou k-mers. Entretanto, estas sequências podem conter bases mal identificadas durante o processo de sequenciamento, de nominadas erros de base calling.

As contigs podem ser usadas para a construção de sequências maiores, denominadas scaffolds. Uma scaffold é resultado da conexão de duas ou mais contigs a partir de uma determinadas evidência (ex: informações de leituras paired-end, posicionamento a partir de um genoma de referência, etc). No caso das leituras paired-end, é possível gerar scaffolds através do alinhamento das leituras contra as contigs e identificação dos casos onde as mates (ou “leituras irmãs”) alinharam com contigs diferentes.

Exemplos de ferramentas para montagem de genomas incluem são: Velvet, Ray, CLC Genome Workbench, MIRA, Newbler, SGA e A5, dentre muitas outras. Destas, a ferramenta A5 é a mais utilizada para dados gerados pela plataforma Illumina MiSeq.

Após um genoma ser montado é necessário se realizar a identificação de seus elementos funcionais (features). Esta etapa é denominada anotação (annotation), e engloba a identificação das regiões de DNA codificante (coding DNA sequences, ou CDSs), regiões de genes de RNAs não-codificantes (ex: tRNAs, rRNAs), anotação funcional das proteínas identificas, dentre outras. Este processo pode ser automatizado com uso de pipelines de anotação, sendo a Prokka a mais utilizada para genomas microbianos.

Por fim, no caso de genomas que serão descritos em artigos científicos, é necessário se disponibilizar as sequências em bancos de dados públicos, como GenBank e o ENA. Para isso, é necessário se formatar o resultado da montagem e anotação em arquivos específicos para a submissão. No caso do GenBank, o formato de submissão é o .sqn.

Neste tutorial partiremos de um conjunto de dados de sequenciamento de genoma completo de Escherichia coli feito com a plataforma Illumina MiSeq, e realizaremos a etapas etapas de montagem, anotação e preparação de arquivo de submissão para o GenBank 😃. Serão utilizados os seguintes programas:

• SRA Toolkit: para baixar os dados de sequenciamento a partir do NCBI SRA.
• A5: montagem do genoma.
• QUAST: avaliação da montagem.
• Prokka: anotação.
• tbl2asn: formatação do arquivo de submissão para o GenBank.

Estrutura do projeto

Este projeto utilizará a seguinte estrutura de arquivos e diretórios:

project/
    bin/
    scripts/
        install_dependencies.sh
        assembly.sh
        assembly_evaluation.sh
        annotation.sh
        generate_submission_file.sh
    reads/
    assembly/
    assembly_evaluation/
    annotation/
    submission/
    Makefile
    samples

Como nos demais projetos, vamos utilizar um arquivo Makefile para centralizar a execução dos scripts de cada etapa. Este arquivo deverá conter o seguinte código:

Instalando as dependências

Como alguns programas que vamos utilizar dependem de uma série de bibliotecas e programas de terceiros, o processo de instalação manual pode se tornar bastante tedioso. Para facilitar as coisas, criei o script install_dependencies.sh, que automatiza esta etapa.

Este script é executado pela regra install_dependencies do Makefile. Para fazer a instalação, utilize o seguinte comando:

$ sudo make install_dependencies

Obtendo o dataset

Neste projeto vamos utilizar dados que serão baixados do NCBI Short Read Archive (SRA), um repositório público de datasets de NGS. Para isso, vamos criar uma arquivo samples na pasta do projeto, contendo o código de acesso da amostra que vamos analisar:

ERR1748534

Agora, vamos criar o script download_data.sh, que é responsável pelo download dos dados. Este script chama o programa fastq-dump do pacote sra-toolkit para cada amostra (neste acaso, apenas uma) listada no arquivo samples.

Neste caso estamos utilizando também a opção --split-3, que indica que nos casos de bibliotecas do tipo paired-end e mate-pair devem ser gerados 2 arquivos: um contendo as leituras da “esquerda” (*_1.fastq) e outro contendo as leituras da “direita” (*_2.fastq).

Para executar este script, use o seguinte comando na pasta do projeto:

$ make download_data

Montando o genoma

O process de montagem será realizado pela ferramenta A5, uma pipeline realiza não apenas o processo de montagem de novo propriamente dito, mas também diferentes etapas de controle de qualidade, tanto de leitura (ex: remoção de adaptadores, trimagem de leituras de baixa qualidade), quanto de contigs e scaffolds (ex: identificação de erros de montagem), como mostrado abaixo:

Etapas de processamento da ferramenta A5

Para executar o A5 vamos utilizar o seguinte código no script assembly.sh:

Note que podemos indicar com os argumentos --begin e --end em qual etapa da pipeline queremos iniciar a nossa análise. Isso é útil caso estejamos re-executando o processo, ou caso não queiramos que alguma etapa final de scaffolding seja realizada, por exemplo. Para executar este script, utilize o seguinte comando na pasta do projeto:

$ make assembly

Avaliação da montagem

A partir do genoma montado podemos derivar diferentes métricas que servem para avaliar a qualidade da montagem gerada. Estas métricas são particularmente úteis para verificarmos o grau de fragmentação do genoma, e dependendo do caso pode ser necessário se realizar um nova montagem, com um novo protocolo, ou mesmo um novo sequenciamento. Para realizar esta análise vamos utilizar a ferramenta QUAST, que será executada pelo script assembly_evaluation.sh.

Para executar este script, use o seguinte comando na pasta do projeto:

$ make assembly_evaluation

Após na pasta assembly_evaluation haverá, para cada amostra indicada um arquivo samples , uma pasta contendo um arquivo report.pdf. Este arquivo contem um sumário com as estatísticas da montagem gerada pelo A5. Destas, algumas das mais relevantes são o número de contigs (neste caso, 140), o tamanho total da montagem (4.638.180 bp) e o valor de N50 (87,306).

O número de contigs está dentro do esperado para um genoma microbiano sequenciado com esta plataforma. O tamanho da montagem está próximo do esperado para um genoma de E. coli, que geralmente possui em torno de 4,7 Mb. Já o valor de N50 indica que 50% das bases no nosso genomas estão presentes em contigs com pelo menos 87.306. Este valor é pouco informativo visto que temos apenas uma montagem, mas é bastante útil quando precisamos escolher a melhor montagem para um organismos após executar diferentes ferramentas e protocolos.

Anotando o genoma

Agora que temos o genoma montado precisamos realizar a identificação dos seus genes e outros elementos funcionais. Esta etapa é realizada pelo script annotation.sh , que é chamado pela regra annotation do Makefile.

Na pasta do projeto, execute o seguinte comando:

$ make annotation

Agora, se dermos um ls na pasta annotation/ERR1748534 veremos que 12 arquivos foram gerados:

$ ls annotation/ERR1748534
annotation.err  annotation.fna  annotation.gff  annotation.tbl
annotation.faa  annotation.fsa  annotation.log  annotation.tsv
annotation.ffn  annotation.gbk  annotation.sqn  annotation.txt

• annotation.err: Indica que há erros na identificação de algumas proteínas, o que exige uma curadoria manual.
• annotation.faa: Sequências de proteínas codificadas pelas das regiões de CDS.
• annotation.ffn: Sequências de DNA das regiões de CDS.
• annotation.fna: Sequências de DNA das scaffolds / contigs usadas como entrada com ID ajustado.
• annotation.fsa: Sequências de DNA das scaffolds / contigs usadas como entrada com ID ajustado e header apropriado para submissão para o GenBank. Usado para a criação do arquivo .sqn.
• annotation.gbk: Sequências + anotação em formato GenBank.
• annotation.gff: Sequências + anotação em formato General Feature Format.
• annotation.log: Log do processo de anotação.
• annotation.sqn: Sequências + anotação em formato de Sequin, usado para submissão para o GenBank. Entretanto, este arquivo não possui todas as informações exigidas para a submissão, sendo necessário gerar um novo (próximo passo).
• annotation.tbl: Anotação em formato feature table. Usado para a criação do arquivo .sqn.
• annotation.tsv: Anotação em formato tab-separated values. Usado para a criação do arquivo .sqn.
• annotation.txt: Sumário da anotação.

Para termos uma visão geral da anotação, vamos ver o conteúdo do arquivo annotation/ERR1748534/annotation.txt:

$ cat annotation/ERR1748534/annotation.txt
organism: Genus species strain 
contigs: 140
bases: 4638180
tRNA: 86
tmRNA: 1
repeat_region: 2
CDS: 4299
rRNA: 13
gene: 4399

Criando o arquivo de submissão para o GenBank

Para prepararmos o arquivo de submissão para o GenBank precisamos de 3 arquivos:

• sequência do genoma em formato FASTA (.fsa);
• Anotação do genoma em formato feature table (.tbl);
• Template de submissão com dados do depositante, autoria e instituição responsável (.sbt).

O arquivo com a sequência do genoma foi gerado no processo de montagem pelo A5, e o arquivo feature table é gerado pelo Prokka. Agora, precisamos criar um arquivo de template para a submissão. Para isso, utilizamos este formulário disponibilizado pelo NCBI. Para cada amostra que temos no arquivo samples vamos criar um arquivo {amostra}.sbt na pasta submission/.

Depois de criar os arquivos, crie o script submission.sh, que executará o programa tbl2asn para cada amostra. Este programa combina as informações dos três arquivos em um único com o formato .sqn (“sequin”), que é utilizado como padrão pelo GenBank na submissão de genes e genomas.

Agora, basta executar este script com o seguinte comando na pasta do projeto:

$ make submission

A partir dos arquivos .sqn é possível fazer a submissão do genoma para o GenBank utilizando o Submission Portal.

obs: não é permitida a submissão de dados obtidos de “terceiros”, sobretudo dados derivados de repositórios públicos. Só podemos submeter novos genomas para o GenBank caso sejamos os pesquisadores responsáveis pelo seu sequenciamento.