Buscar

SEQUENCIAMENTO DE NOVA GERAÇÃO (NGS)

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Você viu 3, do total de 76 páginas

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Você viu 6, do total de 76 páginas

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Você viu 9, do total de 76 páginas

Faça como milhares de estudantes: teste grátis o Passei Direto

Esse e outros conteúdos desbloqueados

16 milhões de materiais de várias disciplinas

Impressão de materiais

Agora você pode testar o

Passei Direto grátis

Prévia do material em texto

Introdução à análise de dados de sequenciadores de nova geração
Versão 2.0.1
Leonardo Varuzza
Abril 2013
2
Sumário
1 Introdução 5
1.1 O que é NGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Como funciona o NGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.1 Preparo da amostra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.2 Amplificação de biblioteca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.3 Sequenciamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.4 Ion Torrent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.5 SOLiD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 Aplicações do NGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.1 Ressequenciamento genômico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.2 Target Sequencing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.3 RNA Seq . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.4 Sequenciamento denovo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3.5 Metagenoma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2 Arquivos de Sequência 15
2.1 Fasta e FastQ Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.1 Fasta format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.2 FastQ format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 SFF File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.1 Converter arquivo SFF para Fasta ou FastQ . . . . . . . . . . . . . . . . . . . . . . 18
2.3 Unmapped BAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 XSQ Format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3 Mapeamento de Sequências 23
3.1 SAM e BAM Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.1 Estrutura do arquivo SAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.2 BAM File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1.3 Samtools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1.4 Picard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Mapeando os reads com o TMAP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.1 Criando o índice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.2 Mapeando os reads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.3 Exemplo: Mapeando os reads de E. coli com o TMAP . . . . . . . . . . . . . . . . 32
3.2.4 Mapeando dados de Long Mate Pair . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3 Bowtie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.1 Utilizando o Bowtie 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.2 Utilizando o Bowtie 2 com o Ion Torrent . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4 Visualizando arquivos BAM com o IGV . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4.1 Importanto o genoma de referência . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3
4 Arquivos de Anotação de Genomas 39
4.1 BED Formats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2 Formatos GFF e GTF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.3 Gerando arquivos de anotação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.3.1 Obtendo anotações do UCSC Browser . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.4 Manipulando arquivos BED . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.4.1 Extrair sequências da regiões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5 Detecção de Variâncias 45
5.1 VCF Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.1.1 Manipulando arquivos VCF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.1.2 Indexando as variantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2 Utilizando o samtools para detectar SNPs . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.2.1 Gerando um arquivo consenso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.3 Utilizando o GATK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3.1 Chamando Variantes no GATK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3.2 Anotando as variantes com o dbSNP . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.4 Utilizando Ion Varriant Caller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.5 Anotando os SNP’s com o snpEff . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6 Montagem denovo 55
6.1 Montando o genoma com o Mira . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.1.1 Montando uma biblioteca de fragmentos utilizando o Mira 3.4 . . . . . . . . . . . . 56
6.1.2 Montando uma biblioteca de mate-pair utilizando o Mira 3.4 . . . . . . . . . . . . 57
6.1.3 Fazendo uma montagem mista com o Mira 3.4 . . . . . . . . . . . . . . . . . . . . 57
6.1.4 Fazendo uma montagem mista com o Mira 3.9 . . . . . . . . . . . . . . . . . . . . 58
6.1.5 Interpretando os resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.1.6 Comparando a montagem com uma referência . . . . . . . . . . . . . . . . . . . . . 60
6.1.7 Visualizando a montagem no Tablet . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7 Apêndices 65
7.1 Ordem dos genótipos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
7.2 Pileup format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
7.3 Samtools VCF file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
7.4 Script para converter os nomes dos cromossomos em um arquivo VCF . . . . . . . . . . . 67
4
Capítulo 1
Introdução
1.1 O que é NGS
Sequenciadores de DNA são equipamentos que leem uma amostra de DNA e geram um arquivo eletrônico
com simbolos que representam a sequência de bases nitrogenadas – A, C, G, T – contidas na amostra. O
primeiro método popular de sequenciamento da DNA foi o de terminação de cadeia de Sanger, publicado
em 1977. Em 1986 foi lançado o primeiro sequenciador automático de DNA, o ABI 370, e em 1998, o
primeiro sequenciador de eletroforese capilar, o ABI 3700. Com a automatização foi possível realizar
grandes projetos de sequenciamento, como o genoma humano, do camundongo e outros. Para realizar
esses projetos foram montados grandes centros com dezenas de máquinas instaladas e ao custo de bilhões
de doláres.
O sequenciamento de nova geração – Next Generation Sequencing, em inglês, ou simplemente NGS
– são processos de sequenciamento de DNA que utilizam metodologias diferentes da de Sanger, com o
objetivo de acelar e baixar o custo do processo de sequênciamento. Apesar de se diferenciarem conside-
ravelmente entre sí todos os sequenciadores de NGS se baseiam no processamento paralelo massivo de
fragmentos de DNA. Enquanto que um sequenciador de eletroforese processa, no máximo, 96 fragmentos
por vez, os sequenciadores de nova geração podem ler até bilhões de fragmentosao mesmo tempo.
Na figura 1.1 vemos a evolução do custo por megabase sequenciada. Nota-se um primeiro decrésimo
em 2004, ano do lançamento do sequenciador 454 da Roche e um decrésimo mais acentuado a partir
de 2006 e 2007, anos em que foram lançados os sequenciadores de nova geração da Illumina e da Life
Technologies. No gráfico, vemos também qual seria a redução de custo, se a tecnologia de sequenciamento
tivesse evoluido segundo a lei de Moore[23]. Vê-se que a evolução do sequenciamento de DNA foi muito
mais acelerada do que dos processadores de computadores. A implicação disso é que os sequenciadores
evoluiram muito mais rápido do que os computadores que analisam os dados gerados, daí a necessidade
computacional para lidar com os dados gerados ter se tornado muito maior do que há 10 anos.
5
yurifroes@outlook.com
Highlight
�����
����
��
���
����
�����
������
����� ����� ����� ����� ����� ����� �����
��
��
��
��
��
��
���
�����������������������������������
�������������
������������
Figura 1.1: Evolução do custo de sequenciamento por megabase. Fonte: http://www.genome.gov/
sequencingcosts/.
Apesar da redução impressionante no custo por megabase, o custo por reação, também chamdado do
custo de apertar o botão start, ainda é bastante alto, da ordem de dezenas de milhares de dólares. Ou
seja, ampliou-se muito a capacidade dos sequenciadores, permitindo até o sequenciamento de mais de um
genoma por corrida, mas sem reduzir muito o custo de operação do equipamento. Um fator limitante na
redução desse custo é o uso de reagentes caros, como bases marcadas por fluoróforos. No final de 2010,
foi lançado o PGM, da Ion Torrent, o primeiro sequenciador a detectar a incorporação dos nucleotídeos
através de um semicondutor, reduzindo, dessa forma, a complexidade do equipamento e o custo da reação
de sequencimento.
1.2 Como funciona o NGS
Cada tecnologia de sequenciamento possui uma estratégia diferente, mas em geral podemos identifi-
car etapas comums entre todos os sequenciadores1: preparo da amostra, amplificação da biblioteca e
sequenciamento.
1.2.1 Preparo da amostra
Primeiro o DNA é fragmentado por um processo químico, mecânico ou enzimática. Cada um desses
fragmentos é chamado de template. Não importa qual o método de fragmentação escolhido, é importante
que ele quebre o DNA de maneira aletória, de forma que todo o genoma seja coberto de maneira o mais
uniforme possível2.
Após a fragmentação, adaptadores, sequências artificiais conhecidas, são incorporados ao template.
Nessa etapa é possível combinar diferentes amostras em uma mesma reação de sequênciamento através
do uso de adaptadores com uma pequena parte , normalmente 5 ou 10 bases, diferente. Depois de ligados
esses adaptadores, as amostras são misturadas, amplificadas e sequenciadas juntas. Após isso, no processo
de sequênciamento, essa parte do adaptador é lida e as amostras são separadas computacionalmente.
Um outro tipode de biblioteca são as de Long Mate-Pair [28], ou LMP. Nesse caso, são gerados
tamplates maiores, da ordem de algumas kilobases. São ligados adaptadores complementares nas pontas
do template e ele é circularizados. Em seguida é feita uma digestão enzimática que gera um único
fragmento com as duas pontas do template separadas adaptador interno (ver figura 1.2). Apesar de ser
um processo mais trabalhoso para o preparo de uma bibliotecas simples de fragmentos, o LMP permite a
detecção de variações estruturais em projetos de ressequenciamento e de geração de scaffolds em projetos
de sequenciamento denovo.
1Nos últimos anos foram propostas tecnologias que tentam eliminar a etapa de amplificação, sequenciando diretamente
uma única molecula de DNA, porém todas as tecnologias propostas sofreram de problemas de baixa acurácia e baixo
throughput.
2Uma excessão à essa regra é o sequenciamento de amplicons, nesse caso o objetivo é sequenciar pequenas regiões que
tem exatamente o tamanho da leitura gerada pelo sequênciador
6
yurifroes@outlook.com
Highlight
yurifroes@outlook.com
Highlight
yurifroes@outlook.com
Highlight
yurifroes@outlook.com
Highlight
Figura 1.2: Processo da amostra de Long Mate Pair.
Um artefato que pode surgir nessa etapa são as duplicações de reads causadas por artefatos de PCR.
Esses artefatos podem gerar distorções na cobertura do genoma e impactar as análises de variações do
genoma ou de expressão do transcriptoma. Por conta disso, as pipelines de análises normalmente possuem
uma etapa em que os reads duplicados são marcados, e consequentemente ignorados nas análises finais.
1.2.2 Amplificação de biblioteca
A amplificação de bibliotecas tem como objetivo gerar em um pequeno espaço físico milhares de cópias
de cada fragmento de DNA produzido na etapa de preparo da amostra. O objetivo dessa amplificação é
aumentar a fonte de sinal luminoso para a maioria dos sequenciadores, e iônico no o caso do Ion Torrent,
que será detectado na etapa de sequenciamento.
O primeiro processo de amplificação desenvolvido para esse propósito foi o PCR de emulsão. Nele
são criados milhões de micro reatores em uma emulsão de óleo. Esses reatores contêm todos os reagentes
necessários para uma reação de PCR e pequenas esferas (também chamdas de beads ou de IonSpheres)
cobertas com a sequência complementar ao adaptador. Esses fragmentos, servem para fixar os clones do
template na esfera e também como primer para a reação de PCR. Ao final do processo a esfera pode
estar em quatro estados:
1. Caso ideal: Um único fragmento foi incorporado à uma única esfera.
2. Beads policlonais: Múltiplos fragmentos foram incorporados a uma esfera.
3. Empty bead: Nenhum fragmento foi incorporado à esfera.
4. Multiple beads: Mais de uma esfera estava presente no micro reator.
O caso de Multiple beads é controlado pelo tamanho do reator, de forma que caiba somente uma esfera
por reator. As empty beads são eliminadas através de uma operação de enriquecimento para beads com
template incorporado. Por fim, as beads policlonais são controladas por meio de um processo estatístico,
que segue uma distriuiçãode Poisson. Basicamente temos muito mais beads do que templates . Por
exemplo, se tivermos 10 vezes mais beads do que templates, espera-se que somente 0,47% das beads
sejam policlonais. Claro que o efeito secundário é que 90% das beads estejam vazias, porém essas beads
são eliminadas por meio do processo de enriquecimento.
O principal fator que afeta a etapa de amplificação é a quantificação do DNA. Se o DNA for subquan-
tificado, ou seja, existe mais DNA na amostra do que o reportado, o resultado vai ser um aumento da
quantidade de beads policlonais. Por outro lado, se o DNA for sobrequantificado, ou seja, se existe menos
DNA do que o reportado, o resultado vai ser uma quantidade muito pequena de beads com fragmentos.
Na tabela 1.1 vemos a porcentagem esperada de beads policlonais e empty beads, se tivermos uma relação
1 para 1 entre beads e fragmentos esperasse que mais de 26% das beads sejam policlonais.
7
Razão Policlonal (%) Empty (%)
10 0.47 90.5
5 1.75 81.9
4 2.65 77.9
3 4.46 71.7
2 9.02 60.7
1 26.4 36.8
Tabela 1.1: Relação entre a razão beads/fragmentos, a probabilidade de beads policlonais e a probabili-
dade de empty beads.
1.2.3 Sequenciamento
O sequenciador é um instrumento que executa uma série de reações químicas. Estes geram sinais que são
detectados e determinam a sequência de bases template se está sendo analisado. A seguir, vamos mostra
o processo de sequenciamento de dois intrumentos que têm abordagens completamente diferentes: o Ion
Torrent e o SOLiD.
1.2.4 Ion Torrent
A maioria dos sequenciadores utiliza uma DNA polimerase para gerar a fita complementar ao template ,
bases marcadas por fluoroforos, e câmeras a detecção. O Ion Torrent é diferente pois a detecção é feita
diretamente. A reação de polimerização gera naturalmenteum um H+3, ou seja, um próton, que altera
o pH do meio. Essa alteração do pH é detectada por um transistor ISFET[1] e convertida em um sinal
elétrico (figura 1.3).
Figura 1.3: Reação de incorporação de uma base pela polimerase.
Outro ponto importante para determinar a sequência é sincronizar a polimerase com a detecção,
tanto no Ion Torrent quanto no 454 essa sincronização é feita pelo controle do tipo disponível para a
polimerase. Por exemplo, suponha que o início do fragmento que se deseja sequenciar seja AGT e que o
sequenciador disponiblize uma certa quantidade de dTTP. A polimerase vai fazer o pareamento do A com
o T e o sinal vai ser detectado pelo transitor ISFET. Para continuar a reação, a polimerase necessita de
um dCTP, porém esse reagente não está disponível e, portanto, a reação para e a leitura da incorporação
é feita. Em seguida, ocorre uma lavagem, e a base seguinte é injetada, e assim por diante em uma série
3A reação de polimerização também gerar um fosfato e esse é o caminho de deteção utilizado pelo 454, com a diferença
que a emissão do fosfato não é detectada diretamente, mas indiretamente através da ativação de uma luciferase que gera
luz
8
de fluxos. Podemos ver na figura 1.4 uma representação dos sinais detectados pelo sensor de um único
poço. É essa informação de intensidade de sinal que é convertida depois na sequência de bases.
Figura 1.4: Flowgram da sequência AATCTTCGT...
Uma questão relevante para os sequenciadores que utilizam fluxos de dNTP’s são os homopolíme-
ros, sequencias contínuas de bases iguais como AAAA, CCCCC e etc4. Nesse caso, todas as bases do
homopolímero vão ser incorporadas em um único fluxo. Felizmente, o sensor ISFET tem uma resposta
bastante linear; portanto, se um A tem um sinal x, um AA vai ter um sinal aproximadamente 2x, e
assim por diante. Na prática, é possível detectar com boa acurácia homopolímeros de até 6 bases.
O último elemento importante, e que diferencia os sequenciadores da nova geração em relação à
anterior, é o paralelismo da reação e da detecção. No Ion Torrent esse paralelismo é obtido pelo uso de
chips de silício. Utilizando o processo CMOS, o mesmo utilizado na fabricação de chips de computador
ou sensores de câmeras digitais, são construídos milhões de poços microscópicos um pouco maiores do
que as esferas com fragmentos de DNA, de forma que, em cada poço, tenha somente uma esfera. No
chip estão também os transitores IsFET que fazem a detecção da mudança de pH, ou seja, cada poço
possui o seu próprio “pH-gâmero” para fazer a detecção do sinal[29].
1.2.5 SOLiD
A sigla SOLiD significa Sequencing by Ligation and Detection e descreve bem o processo de sequencia-
mento utilizado pelo instrumento. Ao invés de utilizar uma polimerase e detectar a incorparação de cada
uma das bases, o SOLiD utiliza octâmeros marcados com fluoróforos para identificar a sequência alvo.
As primeiras 5 bases da probe garantem a especificidade da ligação da probe com o template, enquanto
que as 3 útimas são inosinas que anelam de maneira inespecífica. Conectado à última, inosina temos o
fluróforo que gerará o sinal luminoso a ser detectado pelo sequenciador (ver fig 1.5).
n1 n2 n3 n4 n5 x x x
Fluoróforo
Probe
Figura 1.5: Estrutura da probes utilizadas pelo SOLiD
No SOLiD, assim como no Ion Torrent, cada fragmento é amplificado milhares de vezes na superfície
de uma bead5. Essas beads são então depositadas e fixadas em uma lâmina de vidro. É muito importante
ter essa fixação, porque sabemos que o sinal luminoso que será gerado pelo processo de sequenciamento
está vindo da mesma bead (ou seja, da mesma população de clones geradas de um template) por meio
das coordenadas do ponto luminoso na lâmina.
4Notem que microsatélites com mais de uma base na repetição, como ACACAC, não são homopolimeros
5Esses uma novo modelo do SOLiD, o 5500W, que não utiliza beads. A amplificação dos templates é feita diretamente
na lâmina.
9
A reação de sequenciamento ocorre para cada um dos milhares de clones em cada uma das centenas
de milhões de beads depositadas na lâmina. As etapas dessa reação são, de maneira simplificada, as
seguintes:
1. Na etapa de construção de biblioteca, é adicionado um primer em cada extremidade de cada frag-
mento, chamados de P1 e P2. No processo de sequenciamento é adicionado um primer complementar
à P1, chamado de PA. A última base desse primer alinha com a última base de P1 (fig 1.6 A)
2. É adicionado um pool equimolar de probes. Como temos 4 bases diferentes e uma estrutura de 5
bases mais 3 inosinas, temos portanto 4ˆ5 = 1024 combinações diferentes de probes. As probes vão
se anelar ao longo do template. Após o anelamento, é adicionada uma ligase, que vai fixar somente
a probe que estiver ao lado de uma ponta 5′6. Após a fixação pela ligase ocorre uma lavagem e
todas as probes não fixadas são removidas (fig 1.6 B)
3. É feita a leitura do floróforo. As inosinas, as três últimas bases da probe, são removidas junto com
um floróforo, criando assim uma ponta 5′ livre para fazer a ligação da próxima probe (fig 1.6 C)
4. É feita uma nova incorporação de probes e o processo se repete (fig 1.6 D)
Template
primer
Probe
primer
Template
Probe
primer
Template
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
5’
primer
Template
5’
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
B. Primeira probe anela ao template e
�uóroforo é lido.
A. Primer se liga ao template C. As 3 últimas bases da probe (as inosinas) são 
removidas junto com a probe.
D. Uma nova probe é incorporada e o processo
se repete.
Figura 1.6: Algumas etapas do processo de sequenciamento: ligação do primeiro primer.
A incorporação de probes é repetida 5, 7, 10 ou 15 vezes, dependendo do tamanho desejado de leitura.
Terminado esse ciclos o sistema é aquecido e a fita complementar ao template que foi gerada denatura e
é eliminada. É então incorporado um novo primer, chamado de PB , que alinha uma base a esquerda de
PA. Todo o processo de incorporação de probes é repetido, mas sempre com uma base à esquerda (fig
1.7). Na tabela 1.3 vemos a relação entre primer, ciclo e bases lidas. Vemos que a cada ligação de probe
duas bases são lidas e que cada base é lida por duas probes diferentes. Por exemplo, a base s1 é lida
pelo ciclo 1 do PA e pelo ciclo 2 do PB . Essa construção é chamada de Two Bases Encoding ou 2BE7,
e permite que se faça depois a correção de erros que aumenta a acurácia do processo de sequenciamento
e esta é a principal característica do SOLiD. Vale a pena também resaltar que na segunda ligação, a
última base do primer P1 é sequenciada, isso é muito importante para fazer depois a decodificação das
cores para bases.
6Duas probes consecutivas não ligam porque na extremidade 5′ da probe tem o floróforo, que impede a ligação 5′–3′.
7Existe também uma codificação alternativa chamada de Four Bases Enconding, ou 4BE, que é opcionalmente utilizada
para fazer uma segunda correção de erros chamada de Exact Call Chemistry ou ECC.
10
Cor Fluoróforo
0 Azul FAM
1 Verde Cy3
2 Amarelo TXR
3 Vermelho Cy5
Tabela 1.2: Codificação das cores em números.
Pelo processo de sequenciamento, cada probe 2 de cada 5 bases cobertas por ela, para cobrir todas as
bases duas vezes, temos que utilizar 5 probes, que se alinham em posições diferentes de P1 e permitem
que se cubra todo o fragmento8. Veja a figura 1.8.
Primer Ciclo Bases Lidas
PA 1 s1s2
PA 2 s6s7
PA 3 s11s12
. . . . . . . . .
PB 1 p0s1
PB 2 s5s6
PB 3 s10s11
. . . . . . . . .
Tabela 1.3: Leitura de bases pelas probes do SOLiD
O resultado do sequenciamento é codificado em números de acordo com a tabela 1.2. A relação entre
as duas primeiras bases da probe, n1 e n2, e a cor do fluróforo é dada pela tabela 1.4. Essa tabela
tem diversas propriedades interessantes: ela é simétrica e nenhumacor se repete na mesma linha ou na
mesma coluna (como em um jogo de Soduko). Por causa dessas propriedades temos que, se soubermos
a primeira base do par e a cor fica determinada a segunda base. Suponhamos que a primeira base do
par seja um T, se a cor lida pela probe for verde, a segunda base é, portano, um G. Como o primeira
probe do primer PB lê a última base do adaptar P1, que é conhecida, podemos portanto descobrir qual
é a primeira base da leitura. Tendo a primeira base da leitura e a segunda cor, podemos descobrir a
segunda base, e assim por diante. Podemos pensar nas cores como transformações entre bases, e que se
essas transformações forem encadeads, podemos gerar todas as bases da leitura.
Suponha que a última base de P1 seja um T, e que a seguinte sequência de cores tenha sido obtida:
3 1 3 1 0 2 . Se consultarmos a tabela 1.2, temos que T na primeira base com 3 gera um A, portanto
a primeira base da nossa sequência é um A, o que nos gera o seguinte resultado intermediário:
A 1 3 1 0 2
Vendo agora, a combinação de A com 1 gera um C, e o nosso segundo resultado intermediário é:
AC 3 1 0 2
Continuando a aplicar as transformações chegamos aos seguintes resultados:
ACG 1 0 2
ACGT 0 2
ACGTT 2
ACGTTC
Portanto, a sequência T313102 codificada em color space representa a sequência ACGTTC em base
space.
8Pela construção das probes sequenciaríamos também as 5 últimas bases de P1, o que não é interessante. Por isso, após
o primer PB , é adicionado um espaçador que desloca a posição inicial do sequenciamento 5 bases para frente.
11
primer
Template
Probe
primer
Template
5’
primer
Template
5’
Template
primer
Probe
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
s1 s2p-2... ...p-1 p0 s3 s4 s5 s6s7
B. Como no primer anterior a probe se anela,
porém deslocada uma base à esquerda.
A. Um novo Primer se liga ao template uma
posição para dentro do primer.
C. Novamente as inosinas e o �uóroforo são
removidos.
D. E o processo segue, sempre com uma base
deslocada à esquerda.
Figura 1.7: Algumas etapas do processo de sequenciamento: Ligação do segundo primer.
n2
A C G T
n
1
A 0 1 2 3
C 1 0 3 2
G 2 3 0 1
T 3 2 1 0
Tabela 1.4: Codificação de bases para cores
1.3 Aplicações do NGS
O número de aplicações do NGS é ilimitado. Qualquer coisa que possa ser transformada em DNA pode
ser sequenciada utilizando o mesmo protocolo. Por isso, se alguém quiser criar uma nova análise, basta
modificar a etapa de preparo de biblioteca e de análise e teremos uma nova aplicação. Mesmo assim,
existe um grupo de aplicações que é mais utilizado pela comunidade científica.
1.3.1 Ressequenciamento genômico
Projetos de ressequenciamento têm como objetivo descobrir diferenças entre o genoma de referência e o
genoma de interesse. Os projetos de ressequenciamento normalmente buscam encontrar diferenças entre
o genoma de uma pessoa saudável e o genoma de uma pessoa com alguma doença, como o câncer ou uma
doença hereditária. O ressequenciamento também tem aplicações agropecuárias: pode-se utiliza-lo para
entender a diferença entre genomas de raças de animais, cultivares de plantas ou isolados de bactérias.
O resultado de um projeto de ressequenciamento é a lista das variantes detectadas. Tendo em mãos
essa lista, o pesquisador pode compará-la com os genes anotados do organismo para tentar entender a
relação entre genótipo e fenótipo. Um parâmetro importante para o ressequenciamento é a cobertura
média do genoma, ou seja, quantas vezes em média cada base do genoma foi coberta. Suponha um
genoma de 3 Gbp: para ter uma cobertura de 20× é preciso gerar 60 Gbp de dados. Para detectar
variantes germinativas é recomendada uma covertura entre 20× e 30×, porque essas variantes vão ter
uma frequência em torno de 100% para mutações homozigotas e 50% para variações heterozigotas. Já
para detectar variantes somáticas, é necessária uma cobertura maior, porque essas variantes ocorrem
12
P1
PA
PB
PE
PD
PC
5 10 15 20 250
Figura 1.8: Esquema de cobertura do template pelos probes para uma leitura de 25 bp.
com frequências menores do que as variantes germinativas. A capacidade de detecção da variante vai
depender da frequência mínima que se deseja detectar e da acurácia das leituras geradas.
1.3.2 Target Sequencing
O ressequenciamento de um genoma inteiro gera um volume enorme de dados e permite fazer uma iden-
tificação completa das variações no genoma. Porém, muitas vezes, os pesquisadores só estão interessados
nas variações das regiões codificantes, ou até mesmo num subconjunto de genes e, portanto, sequenciar
todo o genoma para depois selecionar a parte que interessa é bastante ineficiente. Uma alternativa ao
ressequenciamento completo é o target sequencing, na qual a amostra é tratada para selecionar somente
as regiões de interesse.
Existem duas principais abordagens para o target sequencing: captura por hibridização, método utili-
zado pelos kits SureSelect da Agilent e TargetSeq da Life, e captura por amplificação como o RainStorm
da Rain Dance e o AmpliSeq da Life.
Recomenda-se uma alta cobertra da região de interesse pois, como estamos restringindo o sequencia-
mento às regiões de interesse, queremos garantir que todas as bases tenham cobertura suficiente para que
possamos fazer a identificação de variações de maneira confiável. Uma recomendação corrente é cobrir a
região de interesse com uma média de cobertura entre 80 e 100 vezes. Mesmo assim, o target sequencing
oferece vantagens de custo e de processamento devido ao volume reduzido de dados, um genoma humano
sequenciado a 20× gera 60Gbps de sequências, enquanto que um exoma sequenciado a 100× gera 5Gbps
de dados.
1.3.3 RNA Seq
RNA Seq, também chamado the Whole Transcriptome Shotgun Sequencing, é o uso de NGS para sequen-
ciar cDNA com a intenção de capturar a informação do transcriptoma de um organimo. Ao contrário de
outras técnicas, como Microarray ou RT-PCR, o RNA Seq não necessita uma lista pre-definida dos gene
que se deseja detectar. A princípio, qualquer transcrito que estiver sendo expresso pode ser detectado se
o experimento tiver cobertura suficiente. Além disso, o RNA Seq permite detectar, além da expressão,
eventos de splicing alternativo e expressão de genes desconhecidos.
Segundo as diretrizes do projeto ENCODE[27], para estudar somente a expressão de transcritos
polyA, uma cobertura de 20 a 25 milhões de reads curtos mapeados é suficiente. Mas para detectar
transcritos raros ou variações e isoformas em uma amostra de mamíferos, é necessário uma cobertura
de 100 a 200 milhões de reads. O projeto ENCODE também recomenda a realização de duas ou mais
replicatas biológicas. Normalmente não é necessário realizar réplicas técnicas.
Um dos problemas em manipular RNA é a sua fragilidade. Por isso é preciso tomar muito cuidado
na manipulação da amostra e também é preciso avaliar a sua qualidade antes de fazer o sequenciamento
(utilizando, por exemplo, o BioAnalyzer da Agilent). Caso a amostra esteja degradada, o resultado é
13
uma taxa muito baixa de reads mapeados. Outro problema é o RNA Ribossomal, ele corresponde a
uma grande quantidade da massa de RNA de uma célula e, se não for removido da amostra no final, a
maioria dos reads será de RNA ribossomal, o que normalmente não é objetivo do experimento. É possível
utilizar um kit de depleção de RNA Ribossomol, como o Ribominus da Invitrogen, ou então fazer um
enriquecimento para RNA com calda polyA, como o poly(A) Purist também da Invitrogen. Caso haja
interesse também em RNAs não codificantes é melhor utilizar somente o Ribominus. Caso contrário, o
poly(A) Purist é mais eficiente (alguns grupos utilizam os dois para garantir a remoção dos ribossomais).
1.3.4 Sequenciamento denovo
Quando não se tem um genoma de referência, é necessário realizara montagem denovo. O processo de
montagem denovo é muito mais trabalhoso do que o mapeamento com referência. É necessário ter uma
cobertura muito maior e os programas de montagem cometem muito mais erros do que os programas
de mapeamento. Apesar disso, é possível obter bons resultados na montagem de genomas bacterianos
com relativo pouco esforço. Para genomas de organismos superiores, é necessário combinar diversas
tecnologias de sequenciamento e utilizar computadores com quantidades massivas de memória RAM
para efetuar o processo de montagem do genoma.
1.3.5 Metagenoma
A metagenômica é o estudo do material genético extraído diretamente do ambiente. Normalmente,
quando se quer estudar o genoma de uma bactéria, é feita uma cultura para garantir que se está sequen-
ciando um único genoma. Porém, a diversidade de micro-organismos presentes no ambiente é muito maior
do que é possível acessar via sequenciamento individual de bactérias. Por isso, o estudo do metagenoma
é importante.
Existem dois tipos de estudo de metagenomas:
1. Estudo de diversidade utilizando o gene ribossomal 16s: Nesse tipo de estudo amplifica-se por PCR
a sequência 16s e se compara o resultado contra um banco de dados de bactérias conhecidas. Com
isso, é possível avaliar e comparar a diversidade de bactérias presentes na amostra.
2. Shotgun Metagenomics: Nesse segundo caso não se faz nenhuma seleção de alvo. Todo o DNA
extraído da amostra é fragmentado e sequenciado. A análise consiste em montar o “metagenoma”
da amostra para tentar identificar, além da diversidade de genomas, novos genes.
A metagenômica é uma área muito ativa de pesquisa e os métodos de análise ainda são muito manuais.
14
Capítulo 2
Arquivos de Sequência
Em princípio, é bastante simples representar uma sequência de DNA em formato texto. Cada base pode
ser presentada por um caracter: A para Adenina, C para Citosina, G para Guanina e T para timina.
O código oficial para representar DNA é mantido pela IUPAC e inclui também códigos para identificar
bases ambíguas, ou seja, os casos em que não se sabe ao certo a base correta, mas se sabe que deve ser
um C ou T ou algo similar. O código completo está na tabela 2.1.
A Adenina
C Citosina
G Guanina
T (ou U) Timina (ou Uracila)
R A ou G
Y C ou T
S G ou C
W A ou T
K G ou T
M A ou C
B C ou G ou T
D A ou G ou T
H A ou C ou T
V A ou C ou G
N qualquer base
. ou - gap
Tabela 2.1: Código IUPAC para representar o DNA
Além da informação da base, os sequenciadores também produzem uma estimativa da probablidade
da base detectada estar correta. Usualmente essa probabilidade é representada pelo “Phred Quality
Score”, assim chamado porque foi utilizado pela primeira vez no software phred[5]. A fórmula para
calcular os valores de phred é:
Q = −10 log10 Perro (2.1)
Onde Perro é a probabilidade da base ter sido identificada de maneira errada1. Na figura 2.1, temos
o gráfico da função que calcula o phred score e na tabela 2.2 temos exemplos de alguns valors de phred.
Nela vemos que, por exemplo, uma base com phred 20 tem 99% de acurácia, ou seja, uma chance em
100 de estar errada.
1O score phred é construído com base em uma série de preditores de qualidade que são calibrados com dados reais para
gerar uma estimativa de probabilidade de erro.
15
Q Chance de erro Acurácia da base
10 1 em 10 90 %
20 1 em 100 99 %
30 1 em 1000 99,9 %
40 1 em 10000 99,99 %
50 1 em 100000 99,999 %
Tabela 2.2: Valores de qualidade phred, probabilidades de erro e acurácia
 5
 10
 15
 20
 25
 30
 35
 40
 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
Ph
re
d 
Sc
or
e
P de erro
Phred Score
Figura 2.1: Gráfico da função que calcula o Phred Score
Tipicamente, os valores de phred score estão entre 1 e 40. Até é possível sequenciar com uma acurácia
maior que 40, porém, a acurácia das etapas anteriores de 99,99%, e a acurácia do processo inteiro não
pode ser maior do que da parte menos acurada, o que é uma consequência da propagação derros.
A maneira mais simples de representar os valores de qualidade em formato texto é uma lista de
inteiros, o que adiciona de 2 a 3 bytes de dados para cada base sequenciada (1 a 2 caracteres para o valor
da qualidade e um caracter como separador entre os valores). Antes do NGS essa forma de representação
era prática porque o volume de dados era pequeno, mas recentemente se tornou necessário criar novos
formatos que economizassem bytes.
Neste capítulo, vamos analisar os formatos de representação de dados de sequenciamento mais usados
e como manipulá-los utilizando ferramentas gratuitas.
2.1 Fasta e FastQ Files
2.1.1 Fasta format
A forma mais universal de representar sequências biológicas é o formato fasta. Este formato foi criado
para o programa de alinhamento de DNA e proteínas FASTA[21]. Por ser um formato texto muito
simples, ele foi sendo adaptado para as mais diversas aplicações. Mas por outro lado muitas aplicações
não seguem o padrão de maneira estrita.
De maneira geral, um arquivo fasta contém uma série de cabeçalhos seguidos pelo texto da sequência.
Cada cabeçalho é definido por uma linha que se inicia com o caracter >. Pelo padrão, linhas que iniciam
com o caracter ; são consideradas comentários, mas algumas aplicações utilizam o caracter # para essa
função. Abaixo vemos um exemplo de arquivo fasta:
16
>HSBGPG Human gene for bone gla protein (BGP) (fragment)
GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT
ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC
ATGAGAGCCCTCACACTCCTCGCCCTATTGGCCCTGGCCGCACTTTGCATCGCTGGCCAGGCAGGTGAGTGCCCC
CCTGGAGCCCAGGAGGGAGGTGTGTGAGCTCAATCCGGACTGTGACGAGTTGGCTGACCACATCGGCTTTCAGGA
GGCCTATCGGCGCTTCTACGGCCCGGTCTAGGGTGTCGCTCTGCTGGCCTGGCCGGCAACCCCAGTTCTGCTCCT
CTCCAGGCACCCTTCTTTCCTCTTCCCCTTGCCCTTGCCCTGACCTCCCAGCCCTATGGATGTGGGGTCCCCATC
ATCCCAGCTGCTCCCAAATAAACTCCAGAAG
>HSGLTH1 Human theta 1-globin gene (fragment)
CCACTGCACTCACCGCACCCGGCCAATTTTTGTGTTTTTAGTAGAGACTAAATACCATATAGTGAACACCTAAGA
CGGGGGGCCTTGGATCCAGGGCGATTCAGAGGGCCCCGGTCGGAGCTGTCGGAGATTGAGCGCGCGCGGTCCCGG
GATCTCCGACGAGGCCCTGGACCCCCGGGCGGCGAAGCTGCGGCGCGGCGCCCCCTGGAGGCCGCGGGACCCCTG
CTTCTTGCCGTGCTCTCTCGAGGTCAGGACGCGAGAGGAAGGCGC
Neste exemplo, vemos duas sequencias: uma identificada com HSBGPG e outra como HSGLTH1.
Depois de cada identificador temos um espaço e em seguida um comentário sobre a sequência. Depois
disso temos o corpo de cada sequência, o padrão manda que a sequência seja quebrada em linha entre 70
e 132 caracteres. No exemplo, está uma sequência de nucleotídeos, mas o arquivo fasta pode ser utilizado
para representar qualquer tipo de sequência biológica.
O arquivo fasta pode vir acompanhado de um arquivo com extensão .qual. Nesse arquivo estão
os valores de qualidade associados com cada base. Ele tem basicamente o mesmo formato do arquivo
.fasta, porém, no lugar da sequência de bases, ele tem os valores de qualidade na escala phred separados
por espaços. Abaixo temos um exemplo de arquivo .fasta e o respectivo arquivo .qual.
FASTA
>RWBG8:4:5
CTCATTGCCCTCAACACAGTGGAGCGAATTCCTTTGGAAAACCTGCAGATCATCAGAGGAAATATGTACT
ACGAAAATTCCTATGCCTTAGCAGTCTTATCTAACTATGA
QUAL
>RWBG8:4:5
23 16 31 29 30 26 31 31 32 25 32 33 33 29 32 32 31 29 29 29 33 29 31 31 33 31
30 22 29 20 26 18 24 26 12 26 19 25 28 29 13 29 24 32 32 30 30 30 34 33 33 32
33 33 33 33 33 33 29 33 33 25 31 31 29 29 24 24 24 24 24 24 29 29 31 31 18 29
25 29 21 25 25 13 13 15 19 26 22 26 27 28 29 25 31 28 33 29 33 33 33 33 33 29
33 33 34 34 30 30
Como o fasta é um formato de texto, é possível utilizar ferramentas padrão do unix para manipulá-lo.
Para contar o número de sequências em um arquivo fasta, usa-se o seguinte comando:
grep -c "^>" <arquivo.fasta>
Para contar o número de bases no arquivo pode-se utilizar este outro comando:
grep -v "^[>;]" <arquivo.fasta> | wc -c
2.1.2 FastQformat
O formato fasta/qual foi amplamente utilizado para representar sequências geradas por sequenciadores
de Sanger. Porém com o aumento da capacidade dos sequenciadores os arquivos neste formato passaram
a ficar muito pesados. Uma maneira de reduzir o volume de dados foi a criação do formato fastQ. Neste
formato as bases e os valores de qualidade são representados no mesmo arquivo. Além disso as qualidades
são codificadas com caracteres ASCII ao invés de uma sequência de números. Essa codificação é mais
eficiente porque cada valor de qualidade precisa em geral de 3 bytes, dois digitos do valor mais o espaço,
enquanto que no fastQ cada qualidade precisa de somente 1 byte. Abaixo temos um exemplo de arquivo
fastQ:
17
yurifroes@outlook.com
Highlight
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!’’*((((***+))%%%++)(%%%%).1***-+*’’))**55CCF>>>>>>CCCCCCC65
Nesse formato, o cabeçalho é definido pelo caracter @, e é seguido pelo identificador da sequência. Na
linha seguinte tem-se a sequência e, ao contrário do formato fasta, que manda quebrar a sequência em
diversas linhas, no fastQ deve-se ter somente uma linha (mas muitas vezes essa regra não é respeitada).
Na linha seguinte tem-se o caracter + que pode ser seguido pelo identificador da sequência novamente, e
na última linha tem-se os valores de qualidade representados em ASCII.
Para converter o valor de qualidade phred para fastQ adiciona-se 33 a ele e procura-se o caracter
ASCII correspondente. Este é o padrão definido pelo Sanger Institute e mais utilizado atualmente.
Porém as primeiras versões da pipeline de análise da Illumina definem o valor de qualidade como phred
mais 64, o que pode causar problemas para quem estiver manipulando esses arquivos. A partir da versão
1.8 da pipeline de análises a Illumina também passou a adotar o padrão do Sanger.
Para converter de fastQ para fasta pode-se utiliar o programa seqtk que vem junto com o samtools:
seqtk fq2fa in.fastq > out.fasta
2.2 SFF File
O Standard Flow Format[24], ou SFF, é um arquivo binário que foi desenvolvido em conjunto pela
454 Life Science, o Whitehead Institute e pelo Sanger Institute para representar sequências geradas
por um sequenciador baseado em fluxos como o 454 ou o Ion Torrent. Nesse arquivo a sequência é
representada em termos de fluxos, ou seja, as bases e o tamanho normalizado do sinal detectado. Por
exemplo, a sequência AATGG é representada com o equivalente a A3 T1 G2, mas, como a relação
entre o sinal e o tamanho do homopolímero não é totalmente exata, o comprimento dos homopolímeros
é representado por números fracionados. Assim, o exemplo anterior pode ser algo como: A2.9 T1.1
G1.5. Tal informação é perdida devido aos arredondamentos executados quando se converte de SFF
para fasta, mas um programa que trabalhe diretamente com SFF pode utilizar essa informação para
resolver ambiguidades. Essa informação em termos de intensidade de sinal para cada fluxo também é
chamada de flow space em oposição à sequência pura que está em base space.
O arquivo SFF contém um cabeçalho com informações da corrida, assim como um cabeçalho por
read. Além de informar o nome de cada read e o tamanho, esse cabeçalho contém as informações de
clipping de adaptador e clipping de qualidade. Tendo o arquivo SFF, é possível extrair as sequências
sem o clipping, o que pode ser útil para algumas aplicações.
2.2.1 Converter arquivo SFF para Fasta ou FastQ
Muitos programas não lidam diretamente com arquivos SFF, por isso é necessário fazer a conversão para
formatos mais antigos, como fasta ou fastQ2. Na tabela 2.3 está a lista de alguns conversores de SFF
gratuitos.
Nome Linguagem URL
sff_extract Python http://bioinf.comav.upv.es/sff_extract/
sff2fastq C https://github.com/indraniel/sff2fastq
Flower Haskell http://biohaskell.org/
Tabela 2.3: Conversores de SFF
O conversor mais fácil de utilizar é o sff_extract, pois ele é somente um script em python. Porém,
ele é o conversor mais lento. O sff2fastq é o mais rápido, porém é menos flexível, pois só gera arquivo
no formato fastQ. Por fim, o Flower é rápido e flexível, mas necessita ter a linguagem de programação
Haskell instalada no sistema.
2O Torrent Browser fornece o resultado da corrida no formato BAM, mas ele também gera as sequências nos formatos
SFF e fastQ por meio de plugins.
18
yurifroes@outlook.com
Highlight
yurifroes@outlook.com
Highlight
yurifroes@outlook.com
Squiggly
yurifroes@outlook.com
Squiggly
yurifroes@outlook.com
Highlight
yurifroes@outlook.com
Highlight
Converter utilizando o sff_extract
Para utilzar o sff_extract você precisa baixar o script. O nome do arquivo será sff_extract_<versão>.
Renomeie para sff_extract e torne o arquivo executável com o comando chmod +x sff_extract. Para
converter de SFF para fasta/qual utilize o seguinte comando:
sff_extract -o <prefixo> <entrada.sff>
Onde <prefixo> é o nome dos arquivos .fasta e .qual que vão ser gerados e <entrada.sff> é o
nome do arquivo sff que se deseja converter. Este comando vai gerar três arquivos:
• <prefixo>.fasta Com as sequências em formato fasta
• <prefixo>.fasta.qual Com os valores de qualidade em formato phred
• <prefixo>.xml Arquivo XML com as informações extras de cada read, em especial, as informações
de clipping de cada read.
Abaixo temos um exemplo de um registro em um arquivo XML gerado:
<?xml version="1.0"?>
<trace_volume>
<trace>
<trace_name>RWBG8:4:5</trace_name>
<clip_quality_right>149</clip_quality_right>
<clip_vector_left>5</clip_vector_left>
<clip_vector_right>114</clip_vector_right>
</trace>
...
</trace_volume>
Todos os valores de clipping são indexados a partir de 1. Para encontrar a primeira base não trimada,
utiliza-se a seguinte expressão:
first_base_position = max(1, clip_quality_left, clip_adapter_left) (2.2)
E para achar a última base não clipada utiliza-se (caso os clip_adapter_right ou clip_quality_right
sejam 0, o valor é substituído pelo número total de bases no read):
last_base_position = min(clip_quality_right, clip_adapter_right) (2.3)
Portanto, o read no exemplo acima possui bases não clipadas entre as posições 5 e 114. Toda essa
informação de clipagem é indicada por soft clipping, ou seja, as bases do adaptador e as de baixa
qualidade são indicadas em minúsculas. Para fezer o hard clipping e gerar um arquivo fasta “limpo”,
utilize o seguinte comando:
sff_extract -c -o <prefixo> <entrada.sff>
Neste caso, serão gerados os mesmos arquivos, porém o arquivo XML não terá as informações de
clipping porque o clipping já foi feito nos arquivos fasta e qual.
Para gerar o arquivo em formato FastQ, utilize:
sff_extract -Q -o <prefixo> <entrada.sff>
ou
sff_extract -c -Q -o <prefixo> <entrada.sff>
Esse comando vai gerar o arquivo <prefixo>.fastq e <prefixo>.xml, e da mesma forma que no
exemplo anterior, se for adicionada a opção -c, o resultado será com hard clipping.
19
Utilizando o seq_crumbs
O sff_extract foi suplantado pelo pacote seq_crumbs. Diferentemente do programa anterior, o seq_crumbs
foi concebido em módulos, que podem ser combinados de maneira muito mais versátil. Além disso, ele
utiliza a biblioteca biopython. Para baixar o programa, utilize esse endereço:
https://github.com/JoseBlanca/seq_crumbs
Já o biopython está disponível em:
http://biopython.org/wiki/Download
Para fazer a conversão de SFF para fastQ utiliza-se:
sff_exrtract -c -o <output> <entrada.sff>
Note que a versão do seq_crumbs só gera arquivos em formato fastq, não tendo mais a opção de
gerar fasta/qual. Como no programa anterior, a opção -c pede para o programa fazer o trimming das
sequências.
Para fazer a conversão de arquivos SFF de mate-pair vemos uma diferença muito maior, pois temos
que combinar o programa sff_extract com o programa split_matepairs:
sff_extract <entrada.sff> |
split_matepairs -l ION_TORRENT
A combinaçãodesses dois comandos é equivalente ao sff_extract -l. Mas o seq_crumbs oferece
mais opções de manipulação de arquivos. É possível também gerar dois arquivos de mates desentrelaça-
dos, ou seja, um arquivo com a primeira tag e outro com a segunda. Para isso utilize:
sff_extract <entrada.sff> |
split_matepairs -l ION_TORRENT |
pair_matcher -p <orphan.fastq> |
deinterleave_pairs -o <out.1.fastq> <out.2.fastq>
O comando pair_matcher remove as tags que estão sem os seus respectivos pares e coloca no
arquivo <orpha.fastq>, já o deinterleave_pairs separa os pares nos arquivos <out.1.fastq> e
<out.2.fastq>.
Utilizando o Flower
O Flower é um utilitário para ler arquivos SFF feito na linguagem de programação Haskell. Por ser
uma linguagem compilada ele é muito mais do que o sff_extract, que foi escrito em Python. Duas
desvantagens é que o runtime do Haskell é menos ubíquo do que do Python e ele não gera o arquivo
XML com as informações de clipagem necessárias para o Mira (ver 6.1).
Tendo instalado o Haskell, é muito fácil instalar o Flower. Basta utilizar o programa Cabal da seguinte
maneira:
cabal install biosff
Para converter um arquivo SFF para fastQ, utiliza-se:
flower -q <entrada.sff> > saida.fastq
Além de converter o SFF, o flower também permite inspecionar o arquivo. Para visualizar o SFF em
formato texto, use:
flower entrada.sff | more
Por fim, o flower gera uma visualização do flow space em um formato tabulado muito útil. Para
gera-la:
flower -F <entrada.sff>
Segue um exemplo de resultado:
20
...
RWBG8:9:36 81 T 1.05 Qual {unQual = 34}
RWBG8:9:36 82 C 1.03 Qual {unQual = 27}
RWBG8:9:36 83 G 0.00
RWBG8:9:36 84 A 5.72 Qual {unQual = 27},Qual {unQual = 27},
Qual {unQual = 27},Qual {unQual = 27},
Qual {unQual = 27},Qual {unQual = 7}
RWBG8:9:36 85 T 1.93 Qual {unQual = 27},Qual {unQual = 21}
RWBG8:9:36 86 C 0.97 Qual {unQual = 26}
RWBG8:9:36 87 G 0.00
RWBG8:9:36 88 A 0.00
...
Note que mesmo os fluxos que não tiveram sinal são registrados no arquivo SFF e também que,
para os homopolímeros, temos uma estimativa de qualidade para cada uma das bases (quebra de linha
adicionada para facilitar a visualização).
2.3 Unmapped BAM
O formato BAM, discutido na seção sobre mapeamento (seção 3.1.2), é um formato criado para re-
presentar o mapeamento dos reads de NGS em um genoma de referência. Ele é um formato binário,
compactado, indexável e amplamento suportado que pode representar tanto reads mapeados como reads
não mapeados. Os campos opcionais do formato BAM permite grande flexibilidade. É possível adicionar
as informações de flow space do Ion e do 454 ou o color space do SOLiD direto no arquivo BAM, fazendo
dele uma alternativa para os formatos nativos.
Para converter de Unmapped BAM para fastQ, pode-se utizar o comando SamToFastq do Picard
(ver seção 3.1.4) da sequinte maneira:
java -jar SamToFastq.jar I=<input.bam> FASTQ=<out.fastq>
2.4 XSQ Format
O XSQ é o formato nativo do SOLiD 5500. Ele é um formato binário que diminiu em até 60% o número
de bytes necessários para representar o conjunto de sequências no formato fasta/qual. Ele utiliza como
base o HDF5, que é um formato criado para armazenar grandes volumes de dados numéricos.
Esse formato permite representar os reads tanto em base quanto em color space assim como reads
únicos e pareados. Sendo um formato binário, ele também permite o acesso aleatório dos reads contidos
nele. Por fim, é possível adicionar metadados sobre a corrida no próprio arquivo de sequência, facilitando
o armazenamento e a organização dos dados.
As ferramentas de conversão de XSQ para csfasta e vice-versa, assim como a documentação do
formato, estão disponíveis em:
http://www.lifetechnologies.com/us/en/home/technical-resources/software-downloads/xsq-software.
html
Porém, esses conversores só estão disponíveis para Linux CentOS ou RedHat 4 ou 5. Para converter
de XSQ para csfasta utiliza-se:
./convertFromXSQ.sh -f -o <diretorio de sáida> <in.xsq>
Será criado um diretório Libraries dentro do <diretório de saída>. Dentro de Libraries, será
criado um diretório para cada tag – F3, R3, etc – e dentro deste um diretório reads com os arquivos
csfasta, qual e fastq (caso o sequenciamento tenha sido feito com ECC). A opção -f pede para o programa
filtrar os reads marcados como filtrados pelo basecaller e portanto é altamente recomendável.
Também é possível converter de csfasta/qual para XSQ com o comando convertToXSQ.sh. De-
pendendo do tipo de biblioteca, é preciso passar diferentes argumentos. Para biblioteca de fragmentos
utiliza-se:
21
yurifroes@outlook.com
Highlight
yurifroes@outlook.com
Squiggly
./convertToXSQ.sh -x <out.xsq> --mode Fragment
--c1 <F3.csfasta> --q1 <F3.QV.qual>
--libraryName <libname> --runStartTime "0000-00-00 00:00:00"
E para biblioteca de Long Mate Pair :
./convertToXSQ.sh -x <out.xsq> --mode LMP
--c1 <F3.csfasta> --q1 <F3.QV.qual>
--c2 <R3.csfasta> --q2 <R3.QV.qual>
--libraryName <libname> --runStartTime "0000-00-00 00:00:00"
--libraryInsertSizeMinimum <inserto min.>
--libraryInsertSizeMaximum <inserto máx.>
Nesse caso temos que informar, além dos campos do arquivo utilizados no caso da biblioteca de
fragmentos, o arquivos do de sequência e qualidade do mate e, o tamanho mínimo do e máximo do
inserto.
22
Capítulo 3
Mapeamento de Sequências
Após um experimento de NGS é quase sempre necessário fazer o mapeamento das leituras geradas em
um genoma de referência. O primeiro programa para fazer esse mapeamento foi o BLAT1, criado por
Jim Kenth para mapear leituras com alta similaridade contra um genoma de referência. Porém, o BLAT
não consegue alinhar regiões com menos de 40bp e também não tem performance suficiente para alinhar
os milhões de reads gerados por NGS. Os primeiros mapeadores open source otimizados para NGS foram
o MAQ[17] e o SOAP[19], ambos mapeadores que utilizam hash tables para acelerar a busca, com a
diferença que o MAQ indexa os reads enquanto o SOAP indexa o genoma de referência. Com o aumento
do throughput, dos sequenciadores a abordagem de indexar os reads se tornou inviável e por isso os
mapeadores mais recentes indexam a referência.
O problema do uso de hash table é o uso de memória. Alguns mapeadores, como o mapreads
do SOLiD, podem dividir a hash em partes e fazer o mapeamento em etapas de forma a permitir o
mapeamento de grandes genomas mesmo com pouca memória. Alternativamente, é possível utilizar
um FM-Index, que indexa o genoma de referência compactado pela Burrows-Wheeler[33, 6] e entre os
programas que implementam esta técnica estão o BWA[15, 16],bowtie[11] e o SOAPv2[20].
Para aumentar a eficiência os algoritmos citados diminuem a sensibilidade, limitando o número de
mismatches permitidos. Como os reads de NGS normalmente possuem uma qualidade maior no início
do read, utiliza-se uma estratégia de seed and extend: na primeira fase pega-se o início do read, entre 25
e 30 bases, e faz-se a busca permitindo poucos erros – como 2 por exemplo. Encontrada uma posição de
ancoramento do read é feita a extensão do alinhamento de modo a maximizar o score de alinhamento2.
Na seção seguinte é explorado o format SAM/BAM, que é o formato universal para representar
alinhamentos de NGS e nas próximas seções será mostrado como utilizar alguns mapeadores.
3.1 SAM e BAM Files
O formato SAM, Sequence Alignment/Mapping, é um formato texto criado para representar o resultado
do alinhamento dos reads de NGS contra um genoma de referência[8, 18]. A versão binária no SAM
é chamado de BAM, de Binary SAM, e contém exatamente as mesmas informações, mas em formato
binário. Além do formato binário ser mais eficiente, o arquivo BAM também é compactado em blocos,
o que permite ao mesmo tempo a redução do uso de armazemento e o carregamento parcial do arquivo,
o que é muito importante para programas de visualização.
3.1.1 Estrutura do arquivo SAM
Asprimeiras linhas do arquivo SAM são o cabeçalho. Estas linhas começam com o caracter @ seguido
por um código de 2 caracteres que identificam o tipo de informação contina na linha. A lista de códigos
está na tabela 3.1. Depois do cabeçalho temos os alinhamentos propriamente ditos, um por linha, com
campos separados por TABs. Na tabela 3.2 está a lista campos de cada alinhamento.
1O BLAST e o FASTA são programas mais antigos, mas o objetivo deles é fazer buscas em bancos de sequencias e não
especificamente mapear reads em um genoma montado.
2O BWA-long utiliza uma abordagem um pouco diferente. Em vez de fazer o seed and extend, ele realiza um alinhamento
Smith-Waterman contra o índice BWT do genoma
23
Código Descrição
@HD Identifica o início do arquivo SAM e a versão dele.
@SQ Identifica dada uma das sequências no arquivo de
referência
@RG Read group: identifica conjuntos de reads dentro do
arquivo. Esse registro é muito importante quando
são combinados os resultados múltiplas corridas em
um único arquivo.
@PG Lista de programas utilizados no arquivo.
@CO Comentários
Tabela 3.1: Headers do arquivo SAM
Col. Nome Descrição
1 QNAME Nome do read
2 FLAG Bits indicando diversas informações sobre o alinhamento
3 RNAME Nome da sequência na referência
4 POS Posição mais à esquerda do read que se alinha na referência
5 MAPQ Qualidade do Mapeamento
6 CIGAR CIGAR string
7 RNEXT Referência do próximo read no par/segmento
8 PNEXT Posição do próximo read no par/segmento
9 TLEN Tamanho do template observado
10 SEQ Sequência do read
11 QUAL Qualidade do read codificada utilizando a convenção do fastQ do Sanger
12 FIELDs Campos opcionais
Tabela 3.2: Alinhamento no arquivo SAM
Abaixo segue uma descrição mais detalha de cada um dos campos:
1. QNAME: Nome do read ou ’*’ para indicar que a informação não está disponível.
2. FLAG: Bits descrevendo algumas propriedades do alinhamento. o layout dos bits esta abaixo e o
significado de cada bit está na tabela 3.3.
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
M A U UM R RM F L S Q D
Por exemplo, um read corretamente mapeado vai ter um valor de flag 2, mas se ele mapear na
fita oposta a flag vai ser 2 + 16 = 18. Caso o read seja o primeiro de um par a flag vai ser
1 + 2 + 16 + 64 = 83.
3. RNAME: Nome da sequência de referência ou um número indicando qual o registro @SQ do header
com a correspondente sequência. Um read não mapeado tem o valor ’*’ neste campo.
4. POS: Posição mais à esquerda do alinhamento do read com a referência (valores começando em
1). Se o valor for 0, o read não está alinhado.
5. MAPQ: Qualidade do mapeamento. Este valor é representado como um inteiro com interpretação
similar ao índice phred de qualidade, ou seja, P = −10 log10Q, onde Q é a probabilidade do read
estar mapeado de maneira errada.
6. CIGAR: Sequência de caracteres que descrevem como o read está mapeado na referência. O
significado de cada caracter é o seguinte:
24
Sim. valor Descrição
M 1 Template com múltiplos segmentos, ou seja, mate-pair ou pair-
end3.
A 2 Cada segmento está corretamente alinhado.
U 4 Segmento não mapeado.
UM 8 Próximo segmento não mapeado.
R 16 O campo SEQ é o reverso complementar do read.
RM 32 O campo SEQ do próximo segmento é o reverso complementar do
read.
F 64 Esse segmento é o primeiro do template.
L 128 Esse segmento é o último do template.
S 256 Alinhamento secundário. Essa flag indica que existe outro alinha-
mento desse read que é considerado o primário pelo mapeador.
Q 512 Esse alinhamento não passou pelo controle de qualidade.
D 1024 Esse read foi considerado uma duplicata de PCR ou ótica.
Tabela 3.3: Significado e valor de cada flag do format SAM.
M Alignment Match (pode ser um match ou um mismatch)
I Inserção na referência
D Deleção na referência
N Intron (este caracter só é utilizado para representar mRNA’s)
S Soft clipping (corta uma base presente em SEQ)
H Hard clipping (corta uma base não presente em SEQ)
P Padding (deleção em uma padded reference)
= Sequence Match
X Sequence Mismatch
7. RNEXT: Sequência na referência do próximo segmento. Se RNEXT for ’*’ então esta informação
não está presente e se for ’=’ RNEXT tem o mesmo valor e RNAME.
8. PNEXT: Posição do próximo segmento, 0 caso esta informação não esteja disponível.
9. TLEN: Tamanho do inserto entre os segmentos. O segmento à esquerda tem sinal positivo e o
segmento à direita tem sinal negativo, o valor 0 indica que essa informação não está presente.
10. SEQ: Sequência do segmento. Pode ser ’*’ se o valor não tiver sido armazenado. Um ’=’ denota
uma base idêntica à da referência.
11. QUAL: Qualidade das bases de SEQ em formato Sanger FastQ.
Após os campos obrigatórios, pode haver campos opcionais. Esses campos têm o formato
TAG : TY PE : V ALUE
, onde tag é uma sequência de 2 caracteres que inidcam o campo, TYPE indica o tipo de dado e VALUE
é o valor do campo, para a documentação completa veja a especificação do formato SAM[8].
Na figura 3.1 vemos alguns alinhamentos e as entradas equivalentes em formato SAM. O primeiro
read, r001, tem as seguintes flags 163(= 1+2+32+128) = A+M+RM +L, portanto está corretamente
mapeado (A), é o segundo membro de um par (M+L) e o seu par mapeia na posição 37 na fita oposta (flag
RM ). O read r002 posui três soft-clipped bases, a coordenada mostrada no arquivo SAM é da primeira
base alinhada. A string CIGAR deste alinhamento contém um P (padding) que corretamente alinha a
sequência inserida. A informação de padding pode ser omitida se o alinhador não suportar alinhamento
múltiplo de sequências. As últimas 6 bases do read r003 mapeiam na posição 9, e as 5 primeiras mapeiam
na posição 29 da fita reversa. O hard clipping (H) indica bases que não estão presentes na referência. A
tag opcional NM indica o número de mismatches no alinhamento. O read r004 alinha sobre um intron,
fato indicado pelo caracter N.
25
Figura 3.1: Exemplos de alinhamentos e os respectivos registros em formato SAM[18].
3.1.2 BAM File
O arquivo BAM tem exatamente a mesma estrutura do arquivo SAM, porém as informações estão
codificadas em formato binário para economizar espaço em disco e na memória. Além disso, o formato
BAM possui uma compactação em blocos chamada BGZF, que permite o acesso aleatório ao arquivo
mesmo ele estando compactado. Por ter esssa estrutura de blocos, o arquivo BAM pode ser carregado
parcialmente, o que permite a visualização de arquivos BAM muito grandes, sem exceder a capacidade
da memória RAM. Para fazer esse carregamento, o visualizador utiliza um arquivo auxiliar que indexa
os blocos do arquivo BAM em relação às coordenadas genômicas. Esse arquivo tem o mesmo nome do
arquivo BAM, mas com a extensão .bai, de BAM index. Ele também tem que estar no mesmo diretório
que o arquivo BAM.
3.1.3 Samtools
Junto com o formato SAM/BAM foi criada a ferramenta samtools para manipular esses arquivos. Ele
está disponível no site:http://samtools.sourceforge.net/. Para baixar o código fonte do programa,
vá para a URL:
http://sourceforge.net/projects/samtools/files/samtools/0.1.18/
Após baixar o código fonte compile, com o comando make, será gerado um executável chamado
samtools o qual você pode copiar para o diretório que quiser, como por exemplo /usr/local. A
estrutura geral de execução do samtools está abaixo, command é um dos comandos reconhecidos pelo
samtools e <arg>* são as opções para cada um dos comandos.
samtools <command> <arg>*
A seguir, uma lista de receitas úteis para manipular arquivos SAM/BAM:
• Converter entre SAM e BAM: Para fazer a conversão entre a versão binária e texto utiliza-se
o comando view:
– BAM para SAM: samtools view <input.bam> > <output.sam>
– SAM para BAM: samtools view -S -b <input.sam> > <output.bam>
• Visualizar o header de um arquivo BAM: Para extrair somente o header de um arquivo
BAM, use:
samtols view -H <arquivo.bam>• Ordenar e indexar um arquivo BAM: Para visualizar e também para realizar diversas análises
é preciso ordenar os alinhamentos de acordo com a posição genômica. Para isso utiliza-se o seguinte
comando:
26
samtools sort <arquivo.bam> <novo.nome>
Onde <novo.nome> é o prefixo do nome do novo arquivo que vai ser gerado com os reads ordenados.
Esse nome não pode ser o mesmo do arquivo original e não precisa colocar a extensão .bam.
Ela é adicionada automaticamente, e se você colocá-la no final vai ter um arquivo com extensão
.bam.bam, o que não é muito elegante. Caso o arquivo bam não caiba na memória disponível
(que é controlada pela opção -m), será gerado um arquivo temporário com o <novo.nome>.X.bam.
No final do processo esses arquivos vão ser combinados e um arquivo chamado <nome.nome>.bam
é gerado. O valor default de -m é 500000000 (500Mb). Se o computador tiver mais memória
disponível e o arquivo BAM for bastante grande, é interessante aumentar esse valor de forma a
agilizar o processo de ordenação do arquivo bam. Portanto, em um computador com 8Gb de RAM
pode-se utilizar algo como:
samtools sort -m 4000000000 <arquivo.bam> <novo.nome>
Por fim, para indexar o arquivo ordenado, utiliza-se
samtools index <novo.nome.bam>
Que vai gerar o arquivo <novo.nome.bai.
Existe um problema no samtools, o comando sort não altera a tag SO, que indica se o arquivo está
ordenado ou não. Desse modo se o arquivo original tiver SO:unsorted, o arquivo gerado ordena
também vai conter SO:unsorted. Alguns programas estão conscientes desse problemas e lidam
com ele, porém caso seja necessário corrigir essa tag é preciso utilizar o comando reheader do
samtools (veja detalhes sobre esse comando mais abaixo).
• Gerar um relatório sobre o mapeamento:
Para conseguir um relatório dos reads mapeados em arquivo BAM, utiliza-se o comando flagstat
da seguinte maneira:
samtools flagstat <input.bam>
Se for aplicado em um BAM file gerado pelo PGM vai gerar um resultado como este:
99525 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
81284 + 0 mapped (81.67%:-nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (-nan%:-nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (-nan%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
Os valores são classificados como QC-passed reads e QC-failed reads de acordo com o valor da flag
Q(512) do campo FLAGS. Todos os valores estão divididos nessas duas categorias. Esse exemplo
não tem nenhum read marcado com QC-failed e portanto os valores no segundo grupo são todos
zero. A maioria dos valores nesse relatório diz respeito a bibliotecas de reads pareados, e como
essa era uma biblioteca de fragmentos, então o valor mais importante é o número e à porcentagem
de reads mapeados, no caso 81.284, 81,67% do total de 99.525 reads. Abaixo vemos o resultado do
flagstat em um arquivo de reads mate-pair do SOLiD:
27
844482954 + 0 in total (QC-passed reads + QC-failed reads)
225681606 + 0 duplicates
688865129 + 0 mapped (81.57%:nan%)
844482954 + 0 paired in sequencing
422241477 + 0 read1
422241477 + 0 read2
345539940 + 0 properly paired (40.92%:nan%)
533247304 + 0 with itself and mate mapped
155617825 + 0 singletons (18.43%:nan%)
178287164 + 0 with mate mapped to a different chr
90834280 + 0 with mate mapped to a different chr (mapQ>=5)
• Combinar diversos arquivos BAM em um novo arquivo: Para combinar diversos arquivos
bam em um novo arquivo utiliza-se:
samtools merge <out.bam> <input1.bam> <input2.bam> ...
Porém, essa opção vai simplesmente copiar o header do arquivo <input1.bam> e combinar os
alinhamentos de todos os arquivos, mas isso pode fazer com que as tags @SQ, @RG e @PG fiquem
inconsistentes entre o header e os alinhamentos. Por isso, é necessário criar um novo header, com
as informações corretas, e utilizar o comando:
samtools merge -rh <novo.header.txt> <out.bam>
<input1.bam> <input2.bam> ...
Onde <novo.header.txt> é o header gerado pelo usuário. Alternativamente, o usuário pode
utilizar o comando MergeSamFiles do pacote Picard (seção 3.1.4), que atualiza o header ao fazer
o merge.
• Substituir o header de um arquivo BAM: Para substituir o header de um arquivo BAM, crie
um novo header em um arquivo texto e utilize o seguinte comando:
samtools reheader <novo.header.txt> <arquivo.bam> > <novo.arquivo.bam>
Por exemplo, se quiser modificar a tag SO de unsorted para coordinate, utilize os seguintes
comandos:
samtools view -H <arquivo.bam> |
sed ’s/SO:unsorted/SO:coordinate/’ > header.txt
samtools rehader header.txt <arquivo.bam> > <fixed.bam>
3.1.4 Picard
Picard é um conjunto de programas para manipular arquivos SAM/BAM escritos em Java e é um
excelente complemento ao samtools. A página do projeto é http://picard.sourceforge.net/ e um
arquivo .tar.gz, com os pacotes jar já compilados pode ser baixado em:
http://sourceforge.net/projects/picard/files/picard-tools/1.67/
O pacote está estruturado como uma série de arquivos jar, sendo que cada jar é um comando. A
estrutura geral para executar um comando é a seguinte:
java -jar <path de instalação>/<comando>.jar <args*>
Cada ferramenta tem um conjunto próprio de argumentos, porém, alguns argumentos que são comuns
à todas as ferramentas podem ser visualizados na tabela 3.4. Abaixo temos alguns comandos úteis do
Picard, neles o path de instalação está omitido para deixar mais limpa a linha de comando.
28
Argumento Descrição
MAX_RECORD_IN_RAM Número máximo de registros à serem mantidos na
memória. Aumentando esse valor reduze-se o nú-
mero de operações intermediárias de disco. Default:
500.000
TMP_DIR Diretório temporário
CREATE_INDEX Cria arquivo de índice .bai.
Tabela 3.4: Argumentos comuns a todos os comandos do Picard.
1. Ordenar um arquivo BAM: O Picard possui um comando similar ao samtools sort, que é o
SortSam.jar:
java -jar SortSam.jar
I=<entrada.bam> O=<sorted.bam>
SO=coordinate
CREATE_INDEX=true
A grande vantagem do SortSam.jar em relação ao samtools sort é que esse comando atualiza
de maneira apropriada o header do arquivo, trocando a tag SO para coordinate. Para ordenar
arquivos muito grandes é interessante aumentar o valor de MAX_RECORD_IN_RAM. Em uma máquina
com uma boa quantidade de RAM pode-se aumentar esse valor em 10×, como por exemplo:
java -jar SortSam.jar
I=<entrada.bam> O=<sorted.bam>
SO=coordinate
CREATE_INDEX=true
MAX_RECORD_IN_RAM=5000000
2. Fazer o merge de diversos arquivos BAM: O Picard também possui uma alternativa ao
comando samtools merge
java -jar MergeSamFiles.jar
O=<resultado.bam> I=<input1.bam> I=<input2.bam>...
MSD=true
AS=true USE_THREADING=true CREATE_INDEX=true
Assim com o SortSam.jar, o MergeSamFiles.jar trata melhor o header do arquivo gerado do que o
comando equivalente do samtools. Por default ele gera corretamente as tags @RG e @PG no header,
e com a opção MERGE_SEQUENCE_DICTIONARIES, ele também gera as tags @SQ corretas. A opção
AS significa “Assume Sorted” e existe para lidar com o fato do samtools não atualizar o header do
arquivo gerado, e portanto indicar como unsorted um arquivo ordenado. A opção USE_THREADING
faz com que o programa utilize duas threads4.
Quando se deseja combinar diversos arquivos BAM juntos, é um pouco tedioso colocar uma argu-
mento I=<file.bam> para cada arquivo. É possível combinar o MergeSamFiles com o comando
find do unix para combinar todos os arquivos BAM em um mesmo diretório.
java -jar /share/apps/picard/MergeSamFiles.jar
O=<out.bam>
USE_THREADING=true MSD=true
‘find -name \*.bam -printf "I=%p\n"‘
4Linhas de execução do programa. Em um sistema com mais de um processador/core essa opção aumenta a performance
do programa, podendo reduzir o tempo de processamento em 20%
29
3. Identificar Read Duplicados:
Durante o preparo de biblioteca pode ocorrer a duplicaçãode um fragmento, o que causa um viés
artificial da região que originou o fragmento. Por isso, é recomendado que reads suspeitos de serem
duplicatas sejam marcados e ignorados nas análises subsequentes. O MarkDuplicates.jar procura
reads nos quais as posições mais 5’ do alinhamento e a orientação sejam iguais. O read com a
melhor qualidade não é alterado e todos os outros são marcados como duplicados. Para bibliotecas
de mate-pair e pair-end ocorre o mesmo, porém as cordenadas e a orientação dos dois pares de reads
devem coincidir[9]. Uma vantagem do MarkDuplicates.jar é que ele consegue encontrar duplicatas
mesmo se os reads do par não mapearem no mesmo cromossomo.
Para executar o comando utilize:
java -jar MarkDuplicates.jar I=<input.bam>
O=<input.dups.bam>
M=<input.dups_report.txt> AS=true
CREATE_INDEX=true
3.2 Mapeando os reads com o TMAP
O Ion Torrent possui um mapeador otimizado para o padrão de erros dos dados gerados pelo sequen-
ciador. Chamado de TMAP, ele foi criado por Nils Homer. O TMAP é otimizado para os reads de
tamanho variável e também para lidar com erros de homopolímeros gerados pelo Ion Torrent. Para isso,
ele implementa quatro algoritmos de mapeamento que podem ser utilizados individualmente ou combi-
nados. A lista de algoritmos está na tabela 3.5. Na figura 3.2 temos uma comparação da sensitividade
(porcentagem de reads mepados) dos algoritmos do TMAP para o dataset de E. Coli O104H45. Por
comparação foi também colocada a performance bowtie2 e vemos que este tem uma performance um
pouco superior ao map1, porém inferior a todos os outros algoritmos (pelo menos para esse dataset).
Vê-se que o desempenho de map3 e map4 foi bem superior aos outros dois algoritmos, o que gera o
questionamento do porquê deles terem sido incluídos. Na figura 3.3, temos os tempos de processamento
dos mesmos dados da figura anterior. Podemos notar que o map1, apesar de menos eficiente em relação
à porcentagem de reads mapeados, é mais rápido do que map3 e map4.
Para equilibrar sensitividade e performance, temos o comando mapall. Ele aplica cada algoritmo de
mapeamento em estágios definidos pelo usuário. Dessa forma somente os reads que não mapearem em
um estágio são passados para o estágio seguinte. Portanto, uma abordagem eficiênte é colocar o map1
no primeiro estágio, para fazer rapidamente o mapeamento dos reads “fáceis”, e em seguida usar o map3
e o map4 para mapear o resto dos reads.
O código fonte do TMAP está disponível no site github:
https://github.com/iontorrent/TMAP
Para ter uma melhor performance o TMAP pode opcionalmente utilizar o TCMalloc da gperftools
que está disponível em:
http://code.google.com/p/gperftools/
Algoritmo Referência
map1 BWA Short Li and Durbin [15]
map2 BWA Long Li and Durbin [16]
map3 SSAHA Ning [25]
map4 SMEM Li [14]
Tabela 3.5: Algoritmos do TMAP
3.2.1 Criando o índice
O primeiro passo para usar o TMAP é gerar o índice do genoma de referência. Para isso, utiliza-se o
comando index da seguinte forma:
5Disponível na Ion Community
30
bowtie2 map1 map2 map3 map4
60
65
70
75
80
m
a
p 
%
Figura 3.2: Sensitividade dos algortmos do tmap e do bowtie2 (Amostra de E. Coli O104H4)
tmap index -f <referecen.fasta>
Caso haja bases ambíguas na referência, será mostrada uma mensagem de aviso e as bases ambíguas
serão trocadas, dentre as bases possíves, pela de menor valor lexicográfico. Por exemplo, a base Y, que
representa C ou T, será trocada por C, a base N que representa qualquer base será trocada por A, e
assim por diante. A indexação do genoma humano pode demorar 4 horas em um computador atual. O
resultado será diversos arquivos com os índices.
3.2.2 Mapeando os reads
O comando para aplicar um dos algoritmos é:
tmap <mapN> -f <reference> -r <reads>
Onde mapN é map1, map2, map3 ou map4, <reference> é nome de um arquivo fasta que já tenha
sido indexado e <reads> é o arquivo dos reads, que pode estar em formato, fasta, fastQ ou SFF. Como
resultado o tmap gera um arquivo SAM na saída padrão. Caso se utilize um arquivo SFF, é possível
propagar a informação dos flows para o aquivo SAM por meio da opção -Y. Também é possível utilizar
diversas threads com a opção -n <num threads>.
A opção -o X controla o formato de saída do programa, -o 0 é o default e gera uma saída texto em
formato SAM, -o 1 gera um resultado binário em formato BAM, e por fim -o 2 gera um arquivo BAM
não comprimido, uma opção útil se o resultado for tmap vai ser enviado para um outro comando para
mais processamentos, como por exemplo no comando abaixo, que combina o tmap com o samtools sort:
tmap <mapN> -f <reference> -r <reads> -o 2 |
samtools sort -n <memoria> - <output>
Normalmente, os algoritmos individuais não são utilizados, mas sim o comando mapall. Com esse co-
mando você pode combinar os algoritmos em estágios adicionando argumentos no formato stageM mapX*
que define quais algoritmos serão utilizados em cada estágio. É possível definir diversos estágios de modo
que, se um read não mapear em um determinado estágio, ele é passado para o estágio seguinte. Uma
estratégia é colocar o map1 no primeiro estágio, porque ele é mais rápido, e o map3 e o map4 no segundo
estágio. Para isso, utilizamos o seguinte comando:
31
map1 map2 map3 map4
50
10
0
15
0
m
a
p 
tim
e 
(s)
Figura 3.3: Performance dos algoritmos do tmap (Amostra de E. Coli O104H4)
tmap mapall -f <reference> -r <reads>
stage1 map1 stage2 map3 map4
Também é possível controlar as informações da tag @RG através da opção -R. Isso é importante porque
outros programas validam essa tag ao abrir o arquivo BAM. E além disso é bom adicionar metadata
sobre a amostra para facilitar a organização dos arquivos. Algumas tags que podem ser utilizadas são:
ID Identificador único do read group
SM cancer Nome da amostra (sample)
LB hg19 Nome da biblioteca
CN USP Nome do centro de sequencia-
mento
PU PGM/318 A utilização varia entre platafor-
mas. No PGM tem sido usado
para identificar o tipo de chip
utilizado.
PL IONTORRENT Nome da plataforma de sequen-
ciamento (esses valores são tabe-
lados pelo formato SAM)
Combinando todas as informações apresentadas, uma sugestão para rodar o TMAP é a seguinte:
tmap mapall -n <num threads> Y -f <reference> -r <reads.sff>
-R LB:<LB> -R CN:<CN> -R PU:PGM/<chip>
-R ID:<ID> -R SM:<SM> -R PL:IONTORRENT
stage1 map1 stage2 map2 map3 -o 2|
samtools sort - <output>
3.2.3 Exemplo: Mapeando os reads de E. coli com o TMAP
1. Faça o download do arquivo SFF com os reads no site da Ion Community: http://lifetech-it.
hosted.jivesoftware.com/docs/DOC-1518
2. Descompacte o arquivo zip com os reads.
32
3. Faça o download da referência no site do NCBI:http://www.ncbi.nlm.nih.gov/nuccore/NC_
010473.1 em formato fasta e renomeie o arquivo dh10b.fasta
4. Indexe a referência: tmap index -f dh10b.fasta
5. Utilize o TMAP para mapear os reads na referência. Para mapear uma amostra utilize:
tmap mapall -Y -f dh10b.fasta -r 64.sff
-R LB:O104H4 -R CN:uni-munster -R PU:PGM/314
-R ID:64 -R SM:ecoli -R PL:IONTORRENT
stage1 map1 stage2 map2 map3 -o 2 |
samtools sort - 64
Para mapear todas as amostras você pode utilizar o comando for do bash:
for file in *.sff; do
base=‘basename $file .sff‘ ;
tmap mapall -Y -f dh10b.fasta -r $file
-R LB:O104H4 -R CN:uni-munster -R PU:PGM/314
-R ID:$base -R SM:ecoli -R PL:IONTORRENT
stage1 map1 stage2 map2 map3 -o 2|
samtools sort - $base ;
done
6. Marque os reads duplicados: utilize o comando MarkDuplicates.jar do picard para identificar os
reads duplicados:
java -jar MarkDuplicates.jar
I=64.bam O=64.dups.bam M=64.dups.txt AS=true
Ou utilize o for para processar todos os arquivos bam:
for file in *.bam; do
base=‘basename $file .bam‘ ;
java -jar MarkDuplicates.jar
I=$file O=$base.dups.bam M=$base.dups.txt AS=true
done
7. Faça o merge dos arquivos bam utilizando

Outros materiais