Identificando variantes genéticas com NGS

Olá pessoal! Hoje mostarei uma pipeline simples para análise de variantes genéticas a partir de dados de sequenciamento de nova geração (next generation sequencing, NGS).

Variant calling

A análise de variantes genéticas é uma das principais aplicações do NGS hoje em dia, sendo amplamente utilizada na pesquisa básica e aplicada, como o projeto 1000 Genomes. Além disso, esta abordagem já é utilizada também por diversas empresas de diagnóstico molecular não apenas dos Estados Unidos e Europa, mas também no Brasil, como a Mendelics e a Genomika :D

Os principais tipos de variantes que podem ser identificadas com NGS são:

• Single Nucleotide Variants (SNVs): Alterações que resultam na substituição de uma base por outra. Quando ocorrem em uma frequência expressiva em uma determinada população (geralmente ≥ 1% dos individuos), são chamadas Single Nucleotide Polymorphisms (SNPs). O termo SNP é também muitas vezes usado como um sinônimo de SNV.
• Insertions and Deletions (INDELS): Alterações que resultam na inserção ou seleção de uma ou mais bases de uma determinada região do genoma. Quando ocorrem em uma região codificante de DNA (coding DNA sequence, CDS), podem resultar em alterações do quadro de leitura (frameshift).

O processo de identificação de variantes com NGS é geralmente realizado a partir do alinhamento dos dados de sequenciamento contra uma sequência de referência para o organismo de interesse. O resultado deste alinhamento é então processado de modo a se identificar posições do genoma onde há divergência entre o genoma de interesse e o genoma de referência (variant calling).

Neste tutorial vamos identificar variantes em uma cepa de Mycobacterium tuberculosis resistente ao antibiótico Isoniazida, utilizado no tratamento de tuberculose. Para isso, vamos comparar os dados do sequenciamento do genoma desta cepa com o genoma de uma cepa de referência, M. tuberculosis str. H37Rv.

Estrutura do projeto

Para o presente projeto, utilizaremos a seguinte estrutura de arquivos e diretórios:

project/
    scripts/
        download_data.sh
        index_reference.sh
        map_reads.sh
        call_variants.sh
    reference/
        reference.fasta
    reference_indexed/
    reads/
    alignments/
    variants/
    Makefile

Para organizar a execução dos scripts, vamos utilizar um arquivo Makefile, que deverá conter o seguinte código:

Instalando dependências

Neste projeto vamos utilizar os programas:

• SRA-Toolkit: ferramenta para download de dados do NCBI SRA.
• BWA: ferramenta para alinhamento das leituras de sequenciamento contra o genoma de referência (mapeamento).
• Samtools: ferramenta para processamento dos dados de alinhamento.
• Bcftools: ferramenta que realiza o processo de variant calling.

Para instalar estes programas no Ubuntu, utilize o seguinte comando:

$ sudo apt install sra-toolkit bwa samtools bcftools

Obtendo o dataset

O script download_data.sh realizará o download dos arquivos de leitura (em formato .fastq) da cepa de M. tuberculosis resistente a antibióticos a partir do bancos de dados Short Read Archive (SRA) do NCBI com o código de acesso SRR017674.

Para iniciar o download, execute o seguinte comando na pasta raiz do projeto:

$ make download_data

Agora, precisamos do genoma de referência em formato .fasta. Este arquivo pode ser obtido neste link, e deve ser salve na pasta reference.

Indexando a referência

Para que possamos realizar o alinhamento das leituras é necessário primeiro indexar a sequência de referência. Para isso, utilizamos o comando bwa index. O processo de indexação é realizado pelo script index_reference.sh , que é executado pela regra index_reference do Makefile.

Para executar esta etapa, digite o seguinte comando na pasta do projeto:

$ make index_reference

Mapeando as leituras

Após indexar a sequência de referência, podemos realizar o alinhamento das leituras. Para isso, primeiro utilizamos o comando bwa aln para identificar “candidatos” de alinhamentos (regiões com maior probabilidade de serem os sítios de cada leitura), que são armazenados em um arquivo com .sai . Este arquivo é então usado pela comando bwa samse para identificar os alinhamentos mais prováveis, e um arquivo no formato .sam é gerado. Caso estivessemos trabalhando com leituras do tipo paired-end , precisariamos gerar um .sai para cada arquivo e então utilizar o comando sampe para gerar o alinhamento final.

Apesar de conter todos os dados necessários para as etapas posteriores (posição de cada leitura em relação ao genoma de referência), este arquivo precisa ser convertido para um formato com estrutura mais otimizada, a sua “contraparte” binária, o .bam. No nosso caso, esta otimização tem pouco impacto, mas quando trabalhamos com genomas maiores, como o genoma humano, isso acaba tendo um efeito drástico no uso de memória e tempo de processamento.

Para converter o formato .sam para .bam utilizamos o comando samtools view. Depois, é necessário se ordenar as leituras no alinhamento com o comando samtools sort, de modo que estas estejam na ordem em que ocorrem no genoma. Esta etapa permite que o algoritmo de identificação de variantes só precise percorrer o arquivo de alinhamento uma vez. Por fim, o arquivo é indexado com o comando samtools index , e um arquivo .bai é gerado.

Estas etapas são realizadas pelo script map_reads.sh , que é executado pela regra map_reads do Makefile.

Para executar etapa etapa, digite o seguinte comando na pasta do projeto:

$ make map_reads

Identificando variantes

Agora que temos as leituras mapeadas, ordenadas e indexadas, podemos iniciar a análise de variantes. Para isso, primeiro convertemos o formato .bam para o formato de .vcf (variant calling format). Esta etapa é realizada com o comando bcftools mpileup, e nela podemos especificar alguns critérios de filtragem para a nossa análise, com o o uso do argumento -m para indicar qual a “cobertura mínima” de leituras em uma INDEL que deve ser usada como threshold (neste caso, -m 10 é usado). Nesta etapa, geramos um arquiov alignment.vcf.gz.

O arquivo alignment.vcf.gz é então processado pela comando bcftools call para se identificar as regiões com alterações. Como estamos utilizando um genoma procariótico, a opção --ploidy 1 é usada, visto que, por padrão, esta ferramenta considera os genomas como diploides. Este processo gera o arquivo variants.raw.vcf.

Por fim, a ferramenta vcfutils.pl , distribuida com o pacote samtools , é utilizada para se selecionar apenas as variantes com pelo menos 30x de cobertura (argumento -d 30), gerando o arquivo variants.filtered.vcf.

Este processo é realizado pelo script call_Variants.sh.

Para realizar esta etapa, digite o seguinte comando na pasta do projeto:

$ make call_variants

O arquivo variants.filtered.vcf contêm uma tabela com todas as posições do genoma que possuem alterações. Todas as linhas que começam com # são metadados. O campo CHROM (primeira coluna) indica o ID da sequência em que a alteração foi identificada e POS (segunda coluna) indica a posição da sequência onde ocorre a alteração. A base presente na sequência de referência é indicada pela coluna REF (quarta coluna), e a base presente no genoma de interesse é indicada na coluna ALT (quinta coluna). Nos casos de INDELS, o campo INFO (oitava coluna) contêm o campo INDEL; no começo.

A partir das variantes identificas é possível se analisar os seus efeitos com ferramentas como o SnpEff de modo a identificar se as alterações tem efeito biológico (ex: missense, non-sense), e quais genes são afetados. Estas análise podem ser combinadas com a anotação funcional do genoma de referência, permitindo se analisar os processos biológicos que são afetados por estas variantes. É possível também realizar estas análises para diferentes indivíduos, e depois comparar os resultados com informações fenotípicas como é realizado em estudos to tipo Genome Wide Association Studies (GWAS).