Um guia passo a passo para análise de dados de sequenciamento de ADN

Dr. Javier Quilez Oliete, um consultor experiente em bioinformática na Kolabtree, fornece um guia abrangente para análise de dados de sequenciamento de ADN, incluindo ferramentas e software usados para ler dados.

Introdução

Ácido desoxirribonucleico (ADN) é a molécula que transporta a maior parte da informação genética de um organismo. (Em alguns tipos de vírus, a informação genética é transportada pelo ácido ribonucleico (RNA)). Os nucleotídeos (convencionalmente representados pelas letras A, C, G ou T) são as unidades básicas das moléculas de DNA. Conceptualmente, a sequência de ADN é o processo de leitura dos nucleótidos que compõem uma molécula de ADN (por exemplo, “GCAAACCAAT” é uma cadeia de ADN de 10 nucleótidos). As tecnologias atuais de seqüenciamento produzem milhões de leituras de DNA desse tipo em um tempo razoável e a um custo relativamente baixo. Como referência, o custo do sequenciamento de um genoma humano – um genoma é o conjunto completo de moléculas de DNA em um organismo – derrubou a barreira dos 100 dólares e isso pode ser feito em questão de dias. Isto contrasta com a primeira iniciativa de sequenciar o genoma humano, que foi concluída em uma década e teve um custo de cerca de $2,7 bilhões.

Esta capacidade de sequenciar DNA a alto rendimento e baixo custo permitiu o desenvolvimento de um número crescente de métodos e aplicações baseadas em sequenciamento. Por exemplo, o sequenciamento de genomas inteiros ou suas regiões codificadoras de proteínas (duas abordagens conhecidas respectivamente como genoma inteiro e sequenciamento de exomas) em indivíduos doentes e saudáveis pode sugerir alterações de DNA causadoras de doenças. Além disso, o sequenciamento do RNA que é transcrito a partir do DNA – uma técnica conhecida como RNA-sequenciação – é usada para quantificar a actividade genética e como esta muda em diferentes condições (por exemplo, sem tratamento versus tratamento). Por outro lado, os métodos de sequenciamento de captura de conformação cromossômica detectam interações entre moléculas de DNA próximas e assim ajudam a determinar a distribuição espacial dos cromossomos dentro da célula.

Comum a estas e outras aplicações do sequenciamento de DNA é a geração de conjuntos de dados na ordem dos gigabytes e compreendendo milhões de sequências de leitura. Portanto, fazer sentido para os experimentos de seqüenciamento de alta produção (HTS) requer capacidades substanciais de análise de dados. Felizmente, existem ferramentas computacionais e estatísticas dedicadas e fluxos de trabalho de análise relativamente padrão para a maioria dos tipos de dados HTS. Embora algumas das etapas (iniciais) de análise sejam comuns à maioria dos tipos de dados sequenciais, mais análise downstream dependerá do tipo de dados e/ou do objetivo final da análise. Abaixo eu forneço uma cartilha sobre os passos fundamentais na análise de dados HTS e me refiro a ferramentas populares.

Algumas das seções abaixo estão focadas na análise de dados gerados a partir de tecnologias de sequenciamento de leitura curta (principalmente Illumina), uma vez que estas têm historicamente dominado o mercado HTS. Entretanto, novas tecnologias que geram leituras mais longas (por exemplo, Oxford Nanopore Technologies, PacBio) estão ganhando terreno rapidamente. Como o sequenciamento de leitura longa tem algumas particularidades (por exemplo, taxas de erro mais altas), ferramentas específicas estão sendo desenvolvidas para a análise deste tipo de dados.

Controle de qualidade (QC) de leituras brutas

O analista ávido iniciará a análise a partir de arquivos FASTQ; o formato FASTQ tem sido por muito tempo o padrão para armazenar dados de sequenciamento de leitura curta. Em essência, os arquivos FASTQ contêm a seqüência de nucleotídeos e a qualidade por chamada de base para milhões de leituras. Embora o tamanho do arquivo dependa do número real de leituras, os arquivos FASTQ são tipicamente grandes (na ordem de megabytes e gigabytes) e comprimidos. Note que a maioria das ferramentas que usam arquivos FASTQ como entrada podem manuseá-los em formato comprimido, portanto, para economizar espaço em disco, é recomendável não descomprimi-los. Como uma convenção, aqui vou equacionar um arquivo FASTQ a uma amostra sequencial.

FastQC é provavelmente a ferramenta mais popular para realizar o CQ das leituras em bruto. Ele pode ser executado através de uma interface visual ou programática. Enquanto a primeira opção pode ser mais conveniente para usuários que não se sentem confortáveis com o ambiente de linha de comando, a segunda oferece escalabilidade e reprodutibilidade incomparáveis (pense em como pode ser tedioso e sujeito a erros executar manualmente a ferramenta para dezenas de arquivos). De qualquer forma, a principal saída do FastQC é um arquivo HTML com estatísticas resumidas sobre a qualidade geral da sequência em bruto lida a partir de uma determinada amostra. A inspeção de dezenas de relatórios FastQC um por um é tediosa e complica a comparação entre amostras. Portanto, você pode querer usar o MultiQC, que agrega os relatórios HTML do FastQC (assim como de outras ferramentas usadas a jusante, por exemplo, corte do adaptador, alinhamento) em um único relatório.

Informação de QC destina-se a permitir ao usuário julgar se as amostras têm boa qualidade e podem, portanto, ser usadas para os passos subsequentes ou se precisam ser descartadas. Infelizmente, não existe um limiar de consenso baseado na métrica FastQC para classificar as amostras como sendo de boa ou má qualidade. A abordagem que eu utilizo é a seguinte. Espero que todas as amostras que passaram pelo mesmo procedimento (por exemplo, extração de DNA, preparação da biblioteca) tenham estatísticas de qualidade semelhantes e uma maioria de bandeiras de “aprovação”. Se algumas amostras tiverem qualidade inferior à média, ainda as usarei na análise a jusante tendo isto em mente. Por outro lado, se todas as amostras do experimento receberem sistematicamente bandeiras de “aviso” ou “falha” em múltiplas métricas (veja este exemplo), suspeito que algo deu errado no experimento (por exemplo, má qualidade de DNA, preparação da biblioteca, etc.) e recomendo repeti-lo.

Readicionamento de leitura

QC de leituras brutas ajuda a identificar amostras problemáticas, mas não melhora a qualidade real das leituras. Para isso, precisamos de aparar leituras para remover sequências técnicas e fins de baixa qualidade.

Sequências técnicas são restos do procedimento experimental (por exemplo, adaptadores de sequências). Se tais seqüências são adjacentes à seqüência verdadeira da leitura, o alinhamento (veja abaixo) pode mapear leituras para a posição errada no genoma ou diminuir a confiança em um determinado alinhamento. Além das sequências técnicas, também podemos querer remover sequências de origem biológica se estas estiverem muito presentes entre as leituras. Por exemplo, procedimentos de preparação de DNA subótima podem deixar uma alta proporção de RNA ribossômico convertido em DNA (rRNA) na amostra. A menos que este tipo de ácido nucleico seja o alvo do experimento de sequenciamento, manter leituras derivadas do rRNA apenas aumentará a carga computacional das etapas a jusante e pode confundir os resultados. Note que se os níveis de sequências técnicas, rRNA ou outro contaminante forem muito altos, o que provavelmente já terá sido destacado pelo CQ, você pode querer descartar toda a amostra sequenciada.

Na sequência de leitura curta, a sequência de DNA é determinada um nucleotídeo de cada vez (tecnicamente, um nucleotídeo a cada ciclo sequenciado). Em outras palavras, o número de ciclos de seqüenciamento determina o comprimento de leitura. Uma questão conhecida dos métodos de sequenciamento HTS é o decaimento da precisão com que os nucleotídeos são determinados à medida que os ciclos de sequenciamento se acumulam. Isto se reflete em uma diminuição geral da qualidade por chamada de base, especialmente no final da leitura. Como acontece com as seqüências técnicas, tentar alinhar leituras que contêm pontas de baixa qualidade pode levar a uma má colocação ou má qualidade de mapeamento.

Para remover seqüências técnicas/contaminantes e pontas de baixa qualidade, ferramentas de corte de leitura como Trimmomatic e Cutadapt existem e são amplamente utilizadas. Em essência, tais ferramentas irão remover sequências técnicas (internamente disponíveis e/ou fornecidas pelo utilizador) e leituras de corte baseadas na qualidade enquanto maximizam o comprimento de leitura. As leituras que são deixadas muito curtas após o corte são descartadas (leituras excessivamente curtas, por exemplo <36 nucleotídeos, complicam a etapa de alinhamento, já que estas provavelmente mapearão para múltiplos locais no genoma). Você pode querer olhar para a percentagem de leituras que sobrevivem ao corte, já que uma alta taxa de leituras descartadas é provavelmente um sinal de dados de má qualidade.

Finalmente, eu tipicamente reajustei o FastQC nas leituras aparadas para verificar se este passo foi eficaz e melhorou sistematicamente as métricas de QC.

Alinhamento

Com exceções (por exemplo, de novo assembly), o alinhamento (também referido como mapeamento) é tipicamente o próximo passo para a maioria dos tipos de dados e aplicações HTS. O alinhamento de leitura consiste em determinar a posição no genoma do qual a sequência da leitura derivou (tipicamente expressa como cromossoma:start-end). Portanto, neste passo nós exigimos o uso de uma seqüência de referência para alinhar/mapear as leituras em.

A escolha da seqüência de referência será determinada por múltiplos fatores. Para um, a espécie da qual o DNA sequenciado é derivado. Embora o número de espécies com uma sequência de referência de alta qualidade disponível esteja aumentando, este pode ainda não ser o caso para alguns organismos menos estudados. Nesses casos, você pode querer alinhar leituras com uma espécie evolutivamente próxima para a qual um genoma de referência está disponível. Por exemplo, como não há uma sequência de referência para o genoma do coiote, podemos usar a do cão estreitamente relacionado para o alinhamento da leitura. Da mesma forma, podemos ainda querer alinhar nossas leituras com uma espécie intimamente relacionada para a qual existe uma sequência de referência de maior qualidade. Por exemplo, enquanto o genoma do gibão foi publicado, este é dividido em milhares de fragmentos que não recapitulam completamente a organização desse genoma em dezenas de cromossomas; nesse caso, realizar o alinhamento usando a sequência de referência humana pode ser benéfico.

Um outro fator a considerar é a versão do conjunto da sequência de referência, uma vez que novas versões são lançadas à medida que a sequência é atualizada e melhorada. É importante notar que as coordenadas de um determinado alinhamento podem variar entre versões. Por exemplo, múltiplas versões do genoma humano podem ser encontradas no UCSC Genome Browser. Em qualquer espécie, eu sou fortemente a favor de migrar para a mais nova versão assembly uma vez que ela esteja totalmente liberada. Isto pode causar alguns transtornos durante a transição, pois os resultados já existentes serão relativos às versões mais antigas, mas compensa a longo prazo.

Besides, o tipo de dados de seqüenciamento também importa. As leituras geradas a partir dos protocolos DNA-seq, ChIP-seq ou Hi-C serão alinhadas com a sequência de referência do genoma. Por outro lado, como o RNA transcrito do DNA é processado posteriormente em mRNA (ou seja, introns removidos), muitas leituras de RNA-seq falharão no alinhamento com uma seqüência de referência do genoma. Em vez disso, precisamos alinhá-los para transcriptomas seqüências de referência ou usar alinhadores sensíveis à divisão (veja abaixo) ao usar a seqüência do genoma como referência. Relacionado a isto está a escolha da fonte para a anotação da sequência de referência, ou seja, a base de dados com as coordenadas dos genes, transcrições, centrómeros, etc. Eu normalmente uso a anotação GENCODE porque ela combina a anotação abrangente de genes e transcrições de sequências.

Foi desenvolvida uma longa lista de ferramentas de alinhamento de sequências de leitura curta (veja a seção de alinhamento de sequências de leitura curta aqui). Revisá-las está além do escopo deste artigo (detalhes sobre os algoritmos por trás dessas ferramentas podem ser encontrados aqui). Na minha experiência, entre as mais populares estão Bowtie2, BWA, HISAT2, Minimap2, STAR e TopHat. Minha recomendação é que você escolha seu alinhador baseado em fatores chave como o tipo de dados e aplicação HTS, bem como a aceitação pela comunidade, qualidade da documentação e número de usuários. Por exemplo, é necessário alinhar alinhadores como STAR ou Bowtie2 que estejam cientes das junções exon-exon ao mapear o RNA-seq para o genoma.

Comum para a maioria dos mapeadores é a necessidade de indexar a seqüência usada como referência antes que o alinhamento real ocorra. Esta etapa pode ser demorada mas só precisa ser feita uma vez para cada seqüência de referência. A maioria dos mapeadores irá armazenar alinhamentos em arquivos SAM/BAM, que seguem o formato SAM/BAM (arquivos BAM são versões binárias de arquivos SAM). O alinhamento está entre os passos mais demorados e computacionais na análise dos dados de sequenciamento e os arquivos SAM/BAM são pesados (na ordem de gigabytes). Portanto, é importante ter certeza de que você tem os recursos necessários (veja a seção final abaixo) para executar o alinhamento em um tempo razoável e armazenar os resultados. Da mesma forma, devido ao tamanho e formato binário dos arquivos BAM, evite abri-los com editores de texto; ao invés disso use comandos Unix ou ferramentas dedicadas como SAMtools.

Dos alinhamentos

Eu diria que não há um claro passo comum após o alinhamento, já que neste ponto é onde cada tipo de dados e aplicação HTS pode diferir.

Uma análise downstream comum para dados de DNA-seq é a chamada de variante, ou seja, a identificação de posições no genoma que variam em relação à referência do genoma e entre indivíduos. Uma estrutura de análise popular para esta aplicação é o GATK para polimorfismo de nucleotídeo único (SNP) ou pequenas inserções/deleções (indels) (Figura 2). As variantes que compreendem pedaços maiores de DNA (também referidas como variantes estruturais) requerem métodos de chamada dedicados (ver este artigo para uma comparação abrangente). Como com os alinhadores, aconselho selecionar a ferramenta correta considerando fatores chave como o tipo de variantes (SNP, indel ou variantes estruturais), aceitação pela comunidade, qualidade da documentação e número de usuários.

Provavelmente a aplicação mais freqüente do RNA-seq é a quantificação da expressão gênica. Historicamente, as leituras precisavam ser alinhadas à seqüência de referência e então o número de leituras alinhadas a um determinado gene ou transcrição era usado como um proxy para quantificar seus níveis de expressão. Esta abordagem de alinhamento+quantificação é realizada por ferramentas como Cufflinks, RSEM ou featureCounts. Contudo, a abordagem scuh tem sido cada vez mais ultrapassada por novos métodos implementados em software como o Kallisto e o Salmon. Conceptualmente, com tais ferramentas a sequência completa de uma leitura não precisa ser alinhada com a sequência de referência. Ao invés disso, precisamos apenas alinhar nucleotídeos suficientes para ter certeza de que uma leitura se originou de uma determinada transcrição. Simplificando, a abordagem de alinhamento+quantificação é reduzida a um único passo. Esta abordagem é conhecida como pseudo-mapeamento e aumenta muito a velocidade da quantificação da expressão gênica. Por outro lado, tenha em mente que o pseudo-mapeamento não será adequado para aplicações onde o alinhamento completo é necessário (por exemplo, chamada de variante a partir de dados RNA-seq).

Outro exemplo das diferenças nas etapas de análise a jusante e as ferramentas necessárias através da aplicação baseada em sequenciamento é o ChIP-seq. As leituras geradas com tal técnica serão usadas para a chamada de pico, que consiste em detectar regiões no genoma com um excesso significativo de leituras que indicam onde a proteína alvo está ligada. Existem vários picos de chamada e esta publicação os pesquisa. Como exemplo final mencionarei os dados Hi-C, nos quais os alinhamentos são utilizados como input para ferramentas que determinam as matrizes de interação e, a partir destas, as características 3D do genoma. Comentando todos os ensaios baseados em sequências para além do âmbito deste artigo (para uma lista relativamente completa ver este artigo).

Antes de começar…

A parte restante deste artigo toca em aspectos que podem não ser estritamente considerados como passos na análise dos dados HTS e que são largamente ignorados. Em contraste, eu argumento que é capital que você pense sobre as questões colocadas na Tabela 1 antes de começar a analisar dados HTS (ou qualquer tipo de dado de fato), e eu escrevi sobre esses tópicos aqui e aqui.

Tabela 1

Pense sobre isso	Ação proposta
Você tem toda a informação de sua amostra necessária para a análise?	Colher sistematicamente os metadados das experiências
Pode identificar inequivocamente a sua amostra?	Estabelecer um sistema para atribuir a cada amostra um identificador único
Onde estarão os dados e os resultados?	Organização estruturada e hierárquica dos dados
Será capaz de processar várias amostras sem problemas?	Scalabilidade, paralelização, configuração automática e modularidade do código
Poderá você ou qualquer outra pessoa reproduzir os resultados?	Documentar seu código e procedimentos!

Como mencionado acima, os dados brutos do HTS e alguns dos arquivos gerados durante sua análise estão na ordem de gigabytes, portanto não é excepcional que um projeto que inclua dezenas de amostras requeira terabytes de armazenamento. Além disso, alguns passos na análise dos dados HTS são computacionalmente intensivos (por exemplo, alinhamento). Entretanto, a infra-estrutura de armazenamento e computação necessária para a análise dos dados HTS é uma consideração importante e muitas vezes é negligenciada ou não discutida. Como exemplo, como parte de uma análise recente, revisamos dezenas de artigos publicados realizando análise de associação fenomérica (PheWAS). O PheWAS moderno analisa 100-1.000s tanto de variantes genéticas quanto de fenótipos, o que resulta em um importante armazenamento de dados e poder computacional. No entanto, praticamente nenhum dos artigos que revisamos comentou sobre a infra-estrutura necessária para a análise PheWAS. Não surpreendentemente, minha recomendação é que você planeje com antecedência os requisitos de armazenamento e computação que irá enfrentar e os compartilhe com a comunidade.

Need ajuda na análise dos dados de sequenciamento de DNA? Entre em contacto com bioinformáticos freelance e especialistas em genómica em Kolabtree.