Entendo a sopa de letrinhas dos arquivos de sequenciamento genético: FASTQ

 Talvez o título possa parecer um pouco provocativo mas, talvez vocês concordão comigo em breve sobre o formato FASTQ, um dos mais comuns no contexto do sequenciamento massivo em paralelo. 

Os arquivos FASTQ também têm uma estrutura bem definida, como os arquivos FASTA, mas com mais informações.Isso acontece porque eles têm propósitos diferentes: o formato FASTQ surgiu para armazenar dados de leituras de sequenciamento em conjunto com dados de qualidade para cada base lida. As leituras são armazenadas sem uma ordem definida. 

O formato é bastante simples, na verdade, com 4 linhas por read, cada uma com informações diferentes: 

Linha 1: Inicia-se com um @ e em seguida tem o nome da read, além de poder ter alguma descricão

Linha 2: São as bases lidas na read 

Linha 3: Geralmente consiste apenas do caractere +, porém pode repetir o identificador da linha 1

Linha 4: São os dados de qualidade, agora quase sempre calculados como um escore phred, codificados no formato ASCII

Um exemplo de duas reads em um arquivo .FASTQ seria:

  
@201110_A00776_0100_BHKCYNDSXY:1:2550:26521:36307/1
TGGACACTCACATCCTCACCCGTTCCCAATAGGTGCGATGTCCCATCCAGGGATGTCTCAGATGCAGCTGGCAAACCATGGCCCTCATGGCTGAGGACACA
+
?????????????????????????????????????????????????????????????????????????????????????????????????????

@201110_A00776_0100_BHKCYNDSXY:4:2566:19569:4946/1
GATCAGGGTGGTGGCTGCCCCCCAGATCCTGGGGTCCCATCACGGGTGTGTCCTCAGCCATGAGGGCCATGGTTTGCCAGCTGCATCTGAGACATCCCTGG
+
?????????????????????????????????????????????????????????????????????????????????????????????????????

Aqui há algumas questões interessantes que podem ser observadas: o nome das reads tem um padrão, que irá depender do fabricante da plataforma de sequenciamento. Aqui temos uma nomenclatura que é típica da Illumina: o nome do equipamento aparece primeiro (201110_A00776_0100_BHKCYNDSXY), seguido da lane da flowcell (1), além do tile dentro da flowcell onde a read foi lida (2550) e as coordenadas X e Y do cluster dentro do tile (2652:36307), seguido pelo "/1", que representa que estas reads são as primeiras de um par (sequenciada no sentido forward, enquanto as reads 2 são sequenciadas no sentido reverse).

As bases lidas aparecem na segunda linha, e neste caso, não se repete o identificador na linha separadora com o caractere '+'. Talvez vocês estejam se perguntando o que são todos esses pontos de interrogação na quarta linha. São os escores de qualidade das bases lidas, por posição na read. 

Ao ler os escores de qualidade, é importante saber que seguem o formato chamado de Phred, em geral. É o calculo é simples: Q = -10log(p). Aqui, o log está em base 10, e p representa a probabilidade de que a base chamada naquela posição esteja incorreta. Esse número é determinado pelo sequenciador e pelo chamador de bases (que determina qual nucleotídeo foi lido a partir do sinal luminoso, no caso da Illumina).  

É interessante notar que devido à proporção logarítmica, uma mudança 10 no escore para mais representa uma diminuição de 10 vezes na probabilidade de uma chamada incorreta. Assim, um score de qualidade de 10 representa uma chamada incorreta para cada 10 realizadas, com uma acurácia de 90%. O score de 20 já significa uma chamada incorreta para cada 100, com acurácia de 99%, com 30 teríamos 1 incorreta para cada 1000 bases, com acurácia de 99.9%. 

Considerando que o genoma humano é composto por aproximadamente 3.2 bilhões de pares de base, caso queiramos sequenciar todo ele com uma qualidade média de 10 por base, isso significaria cerca de 320 milhões de bases incorretas no sequenciamento, um número bastante alto. Geralmente a qualidade da leitura de 50 é considerada muito boa, e é comum ver nas estatísticas de sequenciamento a porcentagem de bases lidas com qualidade acima de 50 (abreviado como Q50). 

Também é importante notar que o escore de qualidade das bases é diferente de outros escores que veremos mais para frente: de mapeamento e de chamada de variantes. Estes escores levam em consideração a confiabilidade da chamada da base, porém também envolvem outras questões, incluindo a mapeabilidade da região do genoma em que se está alinhando aquela leitura. 

E as interrogações? São os valores de qualidade em formato ASCII. É possível consultar os valores aqui. Há uma certa confusão entre a padronização dos valores, historicamente. O formato definido pelo instituto Sanger é o mais comum hoje em dia (também chamado de fastqsanger). O parâmetro de qualidade vai de 0 a 93 e está representado pelos símbolos ASCII de 33 a 126 (começando com o símbolo ! e terminando com ~ — estes números não fazem sentido, é necessário consultar a tabela). 

A illumina já utilizou diferentes padrões para a codificação da qualidade: inicialmente de 0 a 62 (com o ASCII de 64 a 126). Em algumas versões do padrão da Illumina, a letra B (ASCII 66, representando o valor 2) é utilizada como um indicativo de qualidade baixa de segmento final da leitura. Quando a read termina por bases com qualidade 15 ou menos, a pipeline de chamada substitui o valor de qualidade por 2, que representa que se trata de uma porção da leitura que não é mais confiável e não deve ser utilizada para análise. Quem está familiarizado com a tecnologia da Illumina talvez se lembre que a queda de qualidade nas últimas bases da leitura é uma característica frequente, embora isso tenha melhorado ao longo dos anos.

Por fim, temos um arquivo com milhões de reads não mapeadas, além de seus dados de qualidade. Quando sequenciamos no modo pareado (paired-ends) podemos ter também pares de leituras: Forward e Reverse, ou R1 e R2 que podem ser armazenadas no mesmo arquivo (de modo entrelaçado - uma read 1 seguida pela read 2 correspondente àquele fragmento) ou em 2 arquivos diferentes, que geralmente terão R1 e R2 no final de seu nome. 

O sequenciamento pareado tem a vantagem de permitir obter mais informações sobre o DNA sequenciado. As informações de como as leituras se alinham no genoma (Geralmente é esperado o chamado F1R2 - primeiro forward e depois reverse, porém podemos ver várias variações como F1F2, R1R2 ou R1F2) e do do tamanho do intervalo não sequenciado entre as duas reads (A distância interna do inserto) fornecem evidências sobre possíveis variações estruturais em relação ao genoma de referência: inserções, deleções e inversões. 

Assim fechamos a segunda parte da série. Espero que tenham gostado e em breve teremos uma nova parte.