Alinhamento de
Sequências de Transcritos Gênicos
Daniel Guariz Pinheiro, PhD.
Laboratório de Genética Molecular e Bioinformática
Departamento de Genética
Faculdade de Medicina de Ribeirão Preto
Universidade de São Paulo
Planejamento
• Introdução
• Sequências
– Formatos dos arquivos
• Alinhamento de Sequências
• Algoritmos para Alinhamento de Sequências
• Softwares para Alinhamento de Sequências no Genoma
– Bowtie
– BWA
• Software s auxiliares
– Identificar as sequências
• intersectBED
– Manipular saída de alinhamento (next-generation)
• samtools
– Navegar no genoma
• Integrative Genomics Viewer - IGV
• Exercícios
Introdução
• Mapeamento genômico de transcritos
– Alinhamento das sequências de cDNA obtidas de experimentos (e.g. RNASeq) no genoma referência;
• Identidade aceitável;
– Desafio: qualidade das sequências e configuração dos parâmetros do software de
alinhamento;
• Não ambígua (há uma única identidade);
– Desafio: mapeamentos em regiões repetitivas;
– Identificação das sequências alinhadas em relação a elementos já
mapeados no genoma (e.g. transcritos conhecidos);
• Interseção entre coordenadas de alinhamento;
– Desafio: configuração dos parâmetros para a interseção;
• Identificação de novas regiões com evidências de transcrição;
– Agrupamento de sequências em uma região onde não há transcritos conhecidos;
– Desafio: isolar o que é “ruído”, ou seja, artefatos da técnica;
– Integrar essas informações;
• Relacionar informações de diversas fontes (e.g. expressão gênica, SNPs, regiões
regulatórias, alternative splicing, ...);
• Inspeção visual;
• Utilização de um browser de genoma (e.g. Genome Browser, GBrowse, Integrative
Genomics Viewer (IGV), Gaggle Browser, ... )
Sequências Genéticas
Sequence Read Archive
• SRA (NCBI Sequence Read Archive):
http://www.ncbi.nlm.nih.gov/sra
• ENA (EBI European Nucleotide Archive):
http://www.ebi.ac.uk/ena/
• DRA (DDBJ Sequence Read Archive):
http://trace.ddbj.nig.ac.jp/dra/index_e.shtml
“(…) In mid-September
2010, the
SRA contained >500 billion reads
consisting of 60 trillion base pairs
International
Nucleotide
Sequence
Database
Collaboration
available for download (…) Almost 80% of the
sequencing data are derived from the Illumina
GA platform. The SOLiD™ and Roche/454
platforms account for 15% and 5% of
“We’re growing by about 1 Tb/month.”
submitted base pairs, respectively.(…)”
NCBI’s staff scientist Martin Shumway
(Leinonen R et. al., 2011)
Modelo de Dados
[http://trace.ddbj.nig.ac.jp/dra/documentation_e.shtml]
Formato Fasta
sequence.fa
sequence.qual
>SEQUENCE_1
cagtcagcatactcagtcagtcatgcatgctga
gtcacttgcatgacgtcatgactgcatgactgc
>SEQUENCE_1
1 9 7 15 20 21 16 26 31 37 38 ...
31 13 23 29 31 33 35 30 29 34 ...
Extensões: .fa, .fasta, .fna
Qualidade
O que queremos dizer com qualidade ?
>SEQUENCE_1
1 9 7 15 20 21 16 26 31 37 38 ...
31 13 23 29 31 33 35 30 29 34 ...
Score Perro
Qphred  10log10 (Perror )
10
0.1
20
0.01
30
0.001
Formato fastq
Identificador padrão Illumina
sequence.fastq
@SOLEXA01:1:1:27:1992#0/1
AGTACAAGAGACAGACATTCTTTTTTTTGACACAAG
+SOLEXA01:1:1:27:1992#0/1
\FFFMXPYDDHJSUMVUJLPSNFRXZEDLNLHKHIT
Qualidade codificada como um único caracter da tabela ASCII.
F = 70 (ascii) = 70-64 = 6 (Qphred) = 0,25 (Perror)
Extensões: .fastq
SOLEXA01
the unique instrument name
1
flowcell lane (8 lanes)
1
tile number within the flowcell lane
27
'x'-coordinate of the cluster within the tile
1992
'y'-coordinate of the cluster within the tile
#0
index number for a multiplexed sample (0 for
no indexing)
/1
the member of a pair, /1 or /2 (paired-end or
mate-pair reads only)
P
Qsolexa  10log10 ( error )
1  Perror
Qphred  10log10 (Perror )
Color Space
sequence.csfasta
>9_62_1919_F3
T01231033313313303131110311003330000
T0 = T => T1 = G => G2 = A => ...
sequence.fasta
>9_62_1919_F3
TGATGGCGCATACGCCGTACACCGTGGGCGCCCCC
sff
•
“sff” refere-se a “Standard Flowgram Format”. Os arquivos gerados por uma
corrida de Roche-454.
•
Os arquivos sff contêm:
– um cabeçalho manifesto, no início, o qual descreve informações sobre a corrida;
– valores referentes às intensidades dos sinais para cada base;
•
Formato binário (pirograma) que pode ser convertido para o formato FASTA;
– Programa “sffinfo”.
Sequenciamento em pares
•
Sequenciamento em pares
mate-pair
paired-ends
(Korbel et al. , 2007)
–
–
36 bp
Referência:
>SOLEXA01:1:1:27:1992#0/1
>SOLEXA02:1:1:11:1992#0/1
36 bp
~ ~1928
128bp
bpaa4928
~428bp
bp
paired-ends
mate-pair
>SOLEXA01:1:1:27:1992#0/2
>SOLEXA02:1:1:11:1992#0/2
http://www.illumina.com/technology/mate_pair_sequencing_assay.ilmn
paired-End Sequencing
Mate Pair Library Sequencing for Long Inserts
Sequenciamento em pares (Illumina)
http://www.illumina.com/technology/paired_end_sequencing_assay.ilmn
Introdução
ALINHAMENTO DE SEQUÊNCIAS
Alinhamento de Sequências
Em Bioinformática, alinhamento de sequências é
uma forma de dispor as sequências de DNA,
RNA, ou proteínas para identificar regiões de
similaridade que podem ser consequência de
relacionamentos funcionais, estruturais ou
relações evolutivas entre elas.
Significado Biológico do Alinhamento
de Sequências
• Definição de 3 termos importantes:
– identidade: refere-se à fração de
aminoácidos ou nucleotídeos idênticos entre
pares de sequências após um alinhamento
dessas sequências;
– similaridade: refere-se à fração de
aminoácidos ou nucleotídeos similares (com
propriedades físico-químicas semelhantes –
aminoácidos conservados) entre pares de
sequências após um alinhamento dessas
sequências;
– homologia: representa uma relação
evolutiva entre as sequências;
• Homólogos
– Parálogos;
– Ortólogos;
Há uma referência?
• Resequenciamento
– Existem sequências produzidas a partir de um
genoma/transcriptoma da mesma espécie da amostra
ou de uma espécie relacionada que podem ser usadas
como referências. Alinhamento com a referência.
• Sequenciamento de novo
– Não há sequências que podem ser usadas como
referências. Este tipo de sequenciamento exigirá uma
montagem (assembly) das sequências, utilizando
apenas os dados obtidos desse sequenciamento.
Alinhamento entre as sequencias geradas, que
permitirá a obtenção de um consenso.
Identificação das sequências
• Resequenciamento
– Alinhamento: Conjunto de Sequências X Sequências
Referências (Ex.: Genoma)
>seq1
gcagtcagtcacacatgtca...
>seq2
cgcgcatgcGcgtactctat...
>seq3
tcgagcatcatcagtcgtca...
>seq4
tatgctttatagcgagtcat...
.....
>chrX
atcacacatgtcacatggtcag
ggcatcagtcagtcagtcatgc
gcgcgcatgcCcgtactctatc
tcatgcgtcagtcatgcatgcg
agcagtcatgcatgcatcgcac
tgcatcatacgtcatgcatgaa
.....
Objetivos:
- Eliminar as sequência sem hit
- Eliminar as sequência com hits múltiplos (ambiguous)
- Identificar as sequência com hit único (unambiguous)
Montagem de sequências
• Sequenciamento de novo
– Alinhamentos:
• Conjunto de Sequências X Conjunto de Sequências (alinhamento
pareado)
• Alinhamento Múltiplo de Sequências (MSA)
Consensus :
Seq A
Seq B
Seq C
Seq D
Seq E
Seq F
Seq G
ACAGTACGACAGTACGACCAGTACGATAGCAGTACGATACGACCGA
TCCAGTACGATAGCAGTACGATCAG
GCACAGTACGACCAGTACGATACAGGAAC
CAGGTACGATACGACGGACGGGG
ACAGTACGACAGTACGAAAC
GTACGACCAGTACGATACACT
AACGACAGTACGAAACGGG
TATAGGTACGATACGACGGAC
Introdução
ALGORITMOS PARA ALINHAMENTO
DE SEQUÊNCIAS
Problema básico
• Transformar uma sequência de caracteres em outra:
– Operações:
• inserção
• deleção
• substituição
– Custo de operação:
• Score de substituição
• Penalidade para Gaps (inserção/deleção)
– Qual é a quantidade de operações mínima ?
– Como achar a séries de operações que vai garantir que usamos
a quantidade de operações mínima ?
Exemplo:
Scores:
Match: 2
Mismatch (S): -1
Gap(I): -2
Gap(D): -2
ACGT
||
G-GT
Score (4-2-1): 1
2 matches: 4
1 gap: -2
1 mismatch: -1
Soluções
• Matrix de pontos (dot matrix)
[Goldstein e Gunawardenaa, 2000]
– Informação qualitativa;
Drosophila Dystrobrevin and Mouse ortholog
Soluções
• Matrix de pontos (dot matrix)
– Informação qualitativa;
• Algoritmos de Programação Dinâmica
– Smith-Waterman; Needleman-Wunsch;
• SW é um algoritmo para achar o alinhamento mais
provável com uma estrutura certa;
Alinhamentos de Sequências
•
Alinhamento Global (e.g. Algoritmo de Needleman-Wunsch)
• As sequências envolvidas devem ser alinhadas de um extremo ao outro.
Adequado quando as sequências possuem aproximadamente o mesmo
tamanho.
Seq X : C A T T A G C A G C C T
|
| |
| | |
Seq Y : - A G T A – - A G C - -
•
Alinhamento Local (e.g. Algoritmo de Smith–Waterman)
• Procura-se alinhar apenas as regiões mais similares, independente da
localização relativa de cada região.
Seq X [4,10]: T A G C A G C
| |
| | |
Seq Y [3,7]: T A - - A G C
Alinhamentos (Global/Local) (DNA/Protein)
• FASTA (http://fasta.bioch.virginia.edu/fasta_www2/fasta_list2.shtml)
• EMBOSS Align (http://www.ebi.ac.uk/Tools/emboss/align/)
Matriz de Programação Dinâmica
GG
A
> Score (-2-1): -3
1 gap: -2
1 mismatch: -1
GG
A
> Score(-1-2): -3
1 mismatch: -1
1 gap: -2
GG
> Score(-4-2): -6
2 gaps: -4
1 gap: -2
A
D(i, j) = max
D(i-1, j-1) + s(xi, yj) (diagonal -> match/mismatch)
D(i -1, j) + g (acima -> gap acima)
D(i, j -1) + g (esquerda -> gap esquerda)
Exemplo:
Scores:
Match: 2
Mismatch (S): -1
Gap(I): -2
Gap(D): -2
ACGT
||
G-GT
D(i-1,j-1)
D(i-1,j)
D(i,j-1)
D(i,j)
traceback
Score (4-2-1): 1
2 matches: 4
1 gap: -2
1 mismatch: -1
Solução
• Matrix de pontos (dot matrix)
– Informação qualitativa;
• Algoritmos de Programação Dinâmica
– Smith-Waterman; Needleman-Wunsch;
• SW é um algoritmo para achar o alinhamento mais provável
com uma estrutura certa;
• Por razões de tempo e espaço, não pode ser usado para
alinhamento de sequências de larga escala;
• Utilizações de aproximações (heurísticas);
• Geralmente, quanto mais rápida for a aproximação, mais
distante estará a resposta da solução “correta”;
Desafios
• Eficiência;
– velocidade;
– sensibilidade;
– especificidade;
• Ambiguidade causada por sequências
repetitivas;
• Erros inerentes às técnicas de
sequenciamento.
BLAST
•
•
•
Basic Local Alignment Search Tool
http://blast.ncbi.nlm.nih.gov/
Heurística: dicionário de palavras (hash) – k-word
– Constrói dicionário de palavras de cada k-mer da consulta (score threshold) e as procura nas
referências por uma correspondência exata (match exato);
E-value (S): número de diferentes alinhamentos com scores equivalentes ou
melhores que S que são esperados ocorrer ao acaso em buscas em um banco
de dados aleatório, do mesmo tamanho, com a mesma composição de bases;
QUANTO MENOR... MELHOR!!!
NÃO CONFUNDIR COM P-value (probabilidade)
BLAT
• BLAT—The BLAST-Like Alignment Tool
• http://genome.ucsc.edu/
• Estruturalmente diferente (BLAST)
– Além de outros pontos, o Blat constrói o dicionário de palavras de tamanho k (k-mers)
do banco de dado de sequências referências (database) e faz as buscas nas sequências
as quais se deseja consultar (query);
•
•
•
•
Blat é mais rápido, porém menos sensível;
Possui código especialmente para lidar com intros em alinhamentos RNA/DNA;
Comumente utilizado para localizar uma determinada sequência no genoma ou
determinar a estrutura de exons de um RNA;
Pode ser utilizado para alinhar sequências de Roche/454;
Alinhamento de sequências curtas
• BLAST/BLAT são lentos demais para alinhar
milhões de sequências (Illumina: 35bp100bp/SOLiD: )
• Considerações:
– Não precisamos de um alinhamento sofisticado
como SW;
– Não precisamos de estatísticas com e-value;
– Normalmente, sabemos a quantidade de
mismatches máximas que queremos;
Introdução
NEXT-GENERATION
Introdução à nova geração de
algoritmos de alinhamento
Alinhamentos: construção de índices (referências
ou sequências)
Dependendo da propriedade do índice. Os
algoritmos podem ser classificados em 3
categorias:
– Tabelas hash
– Índices de sufixos/prefixos (Árvores/Arranjos)
– merge-sorting
• slider (Malhis et al., 2009)
Implementações
• Alinhadores para leituras longas (>200bp):
– BLAT, SSAHA2, BWA-SW, gsMapper.
• Alinhadores para leituras curtas:
– Bfast, BioScope, Bowtie, BWA, CLC bio,
CloudBurst, Eland/Eland2, GenomeMapper,
GnuMap, Karma, MAQ, MOM, Mosaik,
MrFAST/MrsFAST, NovoAlign, PASS, PerM, RazerS,
RMAP, SSAHA2, Segemehl, SeqMap, SHRiMP,
Slider/SliderII, SOAP/SOAP2, Srprism, Stampy,
vmatch, ZOOM, ...
Variação na Eficiência
[Lunter et al., 2010]
Introdução à nova geração de
algoritmos de alinhamento
– Tabelas hash
• Paradigma: seed-and-extend
–
–
–
–
–
seed: correspondência (match) exata
tabela hash: índice (chave) é uma string (sequência k-mer)
chave: k-mer (k = 11 – BLAST default)
extensão (Smith Waterman – BLAST)
spaced-seed (match não consecutivo, baseado em um template
111011001110111)
» número de matches do seed é o seu peso (weight)
• Constrói Hash table das leituras e procura na referência
– ELAND (Anthony. J. Cox, 2006, unpublished data)
– MAQ (Li H et al., 2008)
• Constrói Hash table das referências e procura nas leituras
– SOAP (Li R et al., 2008)
Alinhamentos baseados em Hashing
table
• Idéia dos algoritmos de alinhamentos baseados
em hashing tables com spaced-seeds:
read:
agtcgtat
| ||| ||
genoma: acggcacgaggaactcgaatctgacgcatgcagtacta
Se admitirmos 2 mismatches entre a minha sequência e o genoma.
Se separados em 4 fragmentos, vão existir pelo menos 2 fragmentos
sem mismatches, ou seja, com matches exatos !
spaced-seeds
• 6 possibilidades de seeds (templates), com no mínimo
2 fragmentos com match perfeito (ELAND)
read:
agtcgtat
1
agtc----
11110000
2
--tcgt--
00111100
3
----gtat
00001111
4
ag--gt--
11001100
5
--tc--at
00110011
6
ag----at
11000011
Hash Tables
[Chen et al., 2009]
Alinhamento de Sequências com
hashing
• Softwares
– ELAND (Anthony. J. Cox, 2006, unpublished data),
– MAQ (Li H et al., 2008)
– SOAP (Li R et al., 2008)
• Características:
– Para detectar sequências com mais mismatches, precisamos de mais
seeds:
• Mais mismatches => mais tempo
– Algoritmo mais sofisticado para o alinhamento vai requerer mais
tempo:
• Indels/gaps => mais tempo
• Problemas com hashing:
– Memória e tempo
• Precisa de múltiplos processadores e muita memória.
– Necessidade de métodos menos “glutões”
Introdução à nova geração de
algoritmos de alinhamento
– Índices de sufixos/prefixos (Árvores/Arranjos)
• Estrutura de dados que permite a representação de todos os
sufixos/prefixos de um determinada string S. Possibilitando
encontrar as ocorrências de um determinado padrão
(sequência);
• Alinhamento: encontrar os matches exatos no índice,
utilizando representações:
– Árvore de sufixos (Suffix tree)
– Arranjo de sufixos (Suffix array)
» Índice-FM (FM-index) (Ferragina and Manzini, 2000) – BWT
• Alinhadores de sequências:
• Bowtie (Langmead et al., 2009)
• BWA
• BWA-SHORT (Li and Durbin, 2009)
• BWA-SW (Li and Durbin, 2010)
Árvore de sufixos
String indexada: ^GOOGOL$
String procurada: LOL
Mismatches permitidos: 1
Match: GOL
[Li and Durbin, 2009]
Transformação de Burrows–Wheeler
• FM-Index (Ferragina e Manzini, 2000);
– Burrows-Wheeler Transform (BWT)
• Algoritmo usado normalmente em softwares
de compressão (.bzip2);
• Ideia básica: dada uma sequência S de n
símbolos, reordenar os símbolos formando
outra sequência L, que verifica duas condições:
–
–
a probabilidade de um símbolo ser igual ao anterior
é muito elevada;
é possível reconstruir S a partir de L e de mais
alguma informação (primary index);
Transformation T => BWT(T)
Input
^BANANA@
All
Rotations
Sort the
Rows
^BANANA@
@^BANANA
A@^BANAN
NA@^BANA
ANA@^BAN
NANA@^BA
ANANA@^B
BANANA@^
ANANA@^B
ANA@^BAN
A@^BANAN
BANANA@^
NANA@^BA
NA@^BANA
^BANANA@
@^BANANA
Output
BNN^AA@A
Inverse Transformation BWT(T) => T
Input
BNN^AA@A
Add 1
Sort 1
Add 2
Sort 2
B
N
N
^
A
A
@
A
A
A
A
B
N
N
^
@
BA
NA
NA
^B
AN
AN
@^
A@
AN
AN
A@
BA
NA
NA
^B
@^
Add 3
Sort 3
Add 4
Sort 4
BAN
NAN
NA@
^BA
ANA
ANA
@^B
A@^
ANA
ANA
A@^
BAN
NAN
NA@
^BA
@^B
BANA
NANA
NA@^
^BAN
ANAN
ANA@
@^BA
A@^B
ANAN
ANA@
A@^B
BANA
NANA
NA@^
^BAN
@^BA
Add 5
Sort 5
Add 6
Sort 6
BANAN
NANA@
NA@^B
^BANA
ANANA
ANA@^
@^BAN
A@^BA
ANANA
ANA@^
A@^BA
BANAN
NANA@
NA@^B
^BANA
@^BAN
BANANA
NANA@^
NA@^BA
^BANAN
ANANA@
ANA@^B
@^BANA
A@^BAN
ANANA@
ANA@^B
A@^BAN
BANANA
NANA@^
NA@^BA
^BANAN
@^BANA
Add 7
Sort 7
Add 8
Sort 8
BANANA@
NANA@^B
NA@^BAN
^BANANA
ANANA@^
ANA@^BA
@^BANAN
A@^BANA
ANANA@^
ANA@^BA
A@^BANA
BANANA@
NANA@^B
NA@^BAN
^BANANA
@^BANAN
BANANA@^
NANA@^BA
NA@^BANA
^BANANA@
ANANA@^B
ANA@^BAN
@^BANANA
A@^BANAN
ANANA@^B
ANA@^BAN
A@^BANAN
BANANA@^
NANA@^BA
NA@^BANA
^BANANA@
@^BANANA
Output
^BANANA@
FM-Index
Rotações
Transformação de Burrows Wheeler
Sequência
[Li and Durbin, 2009]
Array de sufixos
LF mapping
A ordem de ocorrência de um caracter na última coluna é a mesma ordem de
ocorrência na primeira coluna.
Match exato (aac)
BWT
ponteiros
Suffix Array
Reconstrução da sequência original
(a) T => BWT(T)
(b) BWT(T) => T
(c) Exact Match
Quando dois ponteiros se encontram: match
[Langmead et al., 2009]
Inexact Matching
depth-first-traversal
Bowtie
Backtracks
[Langmead et al., 2009]
BWA-SW
• Heurísticas + Smith-Waterman-like
• Apropriado para sequências longas (> 100 pb)
– Mais acurado e ~10 vezes mais rápido que o BLAT
Sequência Referência
(reference)
Sequência da leitura
(query)
• Percorre as duas estruturas de índices de
sufixos
(para a leitura e para a referência) alinhandoas
utilizando programação dinâmica;
• Utiliza heurísticas para não percorrer
toda a estrutura (Z-best strategy),
identificando seeds de alinhamento com o
maior
score;
representação implícita
A - Árvore de sufixos ‘GOOGOL’;
B – Grafo Direcionado Acíclico de Palavras (DAWG) ‘
GOOGOL’;
• Utiliza o algoritmo de Smith-Waterman para
extender o alinhamento das seeds
identificadas;
Bowtie
• http://bowtie-bio.sourceforge.net
– Burrows-Wheeler;
• Reduz a quantidade de memória e de tempo para
alinhar sequências curtas;
• Podem ser usadas sequências Illumina e SOLiD
– Deficiências:
• Não tem garantia de retornar todos os hits com
mismatches (exceto com opção --best)
• Limite de 3 mismatches (demora mais)
• Reads longos reduz a velocidade
• Não permite alinhamento com gaps
Bowtie Index Builder:
bowtie-build
Usage: bowtie-build [options]* <reference_in> <ebwt_outfile_base>
reference_in
comma-separated list of files with ref sequences
ebwt_outfile_base
write Ebwt data to files with this dir/basename
Options:
-f
reference files are Fasta (default)
-c
reference sequences given on cmd line (as <seq_in>)
-C/--color
build a colorspace index
-a/--noauto
disable automatic -p/--bmax/--dcv memory-fitting
-p/--packed
use packed strings internally; slower, uses less mem
-B
build both letter- and colorspace indexes
--bmax <int>
max bucket sz for blockwise suffix-array builder
--bmaxdivn <int>
max bucket sz as divisor of ref len (default: 4)
--dcv <int>
diff-cover period for blockwise (default: 1024)
--nodc
disable diff-cover (algorithm becomes quadratic)
-r/--noref
don't build .3/.4.ebwt (packed reference) portion
-3/--justref
just build .3/.4.ebwt (packed reference) portion
-o/--offrate <int>
SA is sampled every 2^offRate BWT chars (default: 5)
-t/--ftabchars <int>
# of chars consumed in initial lookup (default: 10)
--ntoa
convert Ns in reference to As
--seed <int>
seed for random number generator
-q/--quiet
verbose output (for debugging)
-h/--help
print detailed description of tool and its options
--usage
print this usage message
--version
print version information and quit
[/data/indexes]$ bowtie-build
$BOWTIE_INDEXES=“/data/indexes”
/data/hg18.fa
hg18
hg18.1.ebwt
hg18.2.ebwt
hg18.3.ebwt
hg18.4.ebwt
hg18.rev.1.ebwt
hg18.rev.2.ebwt
Bowtie Index Inspector:
bowtie-inspect
Usage: bowtie-inspect [options]* <ebwt_base>
<ebwt_base>
ebwt filename minus trailing .1.ebwt/.2.ebwt
By default, prints FASTA records of the indexed nucleotide sequences to
standard out. With -n, just prints names. With -s, just prints a summary of
the index parameters and sequences. With -e, preserves colors if applicable.
Options:
-a/--across <int>
-n/--names
-s/--summary
-e/--ebwt-ref
-v/--verbose
-h/--help
--help
Number of characters across in FASTA output (default: 60)
Print reference sequence names only
Print summary incl. ref names, lengths, index properties
Reconstruct reference from ebwt (slow, preserves colors)
Verbose output (for debugging)
print detailed description of tool and its options
print this usage message
[/data/indexes]$ bowtie-inspect
-s
hg18
Bowtie Aligner: bowtie
Usage:
bowtie [options]* <ebwt> {-1 <m1> -2 <m2> | --12 <r> | <s>} [<hit>]
<m1>
<m2>
<r>
<s>
Comma-separated list of files containing upstream mates (or the
sequences themselves, if -c is set) paired with mates in <m2>
Comma-separated list of files containing downstream mates (or the
sequences themselves if -c is set) paired with mates in <m1>
Comma-separated list of files containing Crossbow-style reads. Can be
a mixture of paired and unpaired. Specify "-" for stdin.
Comma-separated list of files containing unpaired reads, or the
sequences themselves, if -c is set. Specify "-" for stdin.
File to write hits to (default: stdout)
<hit>
Input:
-q
-f
-r
-c
-C
-Q/--quals <file>
--Q1/--Q2 <file>
-s/--skip <int>
-u/--qupto <int>
-5/--trim5 <int>
-3/--trim3 <int>
--phred33-quals
--phred64-quals
--solexa-quals
--solexa1.3-quals
--integer-quals
Alignment:
-v <int>
or
-n/--seedmms <int>
-e/--maqerr <int>
-l/--seedlen <int>
--nomaqround
-I/--minins <int>
-X/--maxins <int>
--fr/--rf/--ff
--nofw/--norc
--maxbts <int>
--pairtries <int>
-y/--tryhard
--chunkmbs <int>
query input files are FASTQ .fq/.fastq (default)
query input files are (multi-)FASTA .fa/.mfa
query input files are raw one-sequence-per-line
query sequences given on cmd line (as <mates>, <singles>)
reads and index are in colorspace
QV file(s) corresponding to CSFASTA inputs; use with -f -C
same as -Q, but for mate files 1 and 2 respectively
skip the first <int> reads/pairs in the input
stop after first <int> reads/pairs (excl. skipped reads)
trim <int> bases from 5' (left) end of reads
trim <int> bases from 3' (right) end of reads
input quals are Phred+33 (default)
input quals are Phred+64 (same as --solexa1.3-quals)
input quals are from GA Pipeline ver. < 1.3
input quals are from GA Pipeline ver. >= 1.3
qualities are given as space-separated integers (not ASCII)
report end-to-end hits w/ <=v mismatches; ignore qualities
max mismatches in seed (can be 0-3, default: -n 2)
max sum of mismatch quals across alignment for -n (def: 70)
seed length for -n (default: 28)
disable Maq-like quality rounding for -n (nearest 10 <= 30)
minimum insert size for paired-end alignment (default: 0)
maximum insert size for paired-end alignment (default: 250)
-1, -2 mates align fw/rev, rev/fw, fw/fw (default: --fr)
do not align to forward/reverse-complement reference strand
max # backtracks for -n 2/3 (default: 125, 800 for --best)
max # attempts to find mate for anchor hit (default: 100)
try hard to find valid alignments, at the expense of speed
max megabytes of RAM for best-first search frames (def: 64)
Reporting:
-k <int>
-a/--all
-m <int>
-M <int>
--best
--strata
Output:
-t/--time
-B/--offbase <int>
--quiet
--refout
--refidx
--al <fname>
--un <fname>
--max <fname>
--suppress <cols>
--fullref
Colorspace:
--snpphred <int>
or
--snpfrac <dec>
--col-cseq
--col-cqual
--col-keepends
SAM:
-S/--sam
--mapq <int>
--sam-nohead
--sam-nosq
--sam-RG <text>
Performance:
-o/--offrate <int>
-p/--threads <int>
--mm
--shmem
Other:
--seed <int>
--verbose
--version
-h/--help
report up to <int> good alignments per read (default: 1)
report all alignments per read (much slower than low -k)
suppress all alignments if > <int> exist (def: no limit)
like -m, but reports 1 random hit (MAPQ=0); requires --best
hits guaranteed best stratum; ties broken by quality
hits in sub-optimal strata aren't reported (requires --best)
print wall-clock time taken by search phases
leftmost ref offset = <int> in bowtie output (default: 0)
print nothing but the alignments
write alignments to files refXXXXX.map, 1 map per reference
refer to ref. seqs by 0-based index rather than name
write aligned reads/pairs to file(s) <fname>
write unaligned reads/pairs to file(s) <fname>
write reads/pairs over -m limit to file(s) <fname>
suppresses given columns (comma-delim'ed) in default output
write entire ref name (default: only up to 1st space)
Phred penalty for SNP when decoding colorspace (def: 30)
approx. fraction of SNP bases (e.g. 0.001); sets --snpphred
print aligned colorspace seqs as colors, not decoded bases
print original colorspace quals, not decoded quals
keep nucleotides at extreme ends of decoded alignment
write hits in SAM format
default mapping quality (MAPQ) to print for SAM alignments
supppress header lines (starting with @) for SAM output
supppress @SQ header lines for SAM output
add <text> (usually "lab=value") to @RG line of SAM header
override offrate of index; must be >= index's offrate
number of alignment threads to launch (default: 1)
use memory-mapped I/O for index; many 'bowtie's can share
use shared mem for index; many 'bowtie's can share
seed for random number generator
verbose output (for debugging)
print version information and quit
print this usage message
[/data]$ bowtie hg18 \
> -c "AGGAATTGCGGGAGGAAAATGGGTAGTTAGCTATTT,AGGGCCCATAGCAACAGATTTCTAGCCCCCTGAAGA"
> --best --strata --tryhard -m 1
Principais parâmetros do Bowtie
Alignment:
-n <int>: número máximo de mismatches na seed [1..3] (2) – Mutuamente exclusivo (-v);
-v <int>: número máximo de mismatches em todo o alinhamento, ignorando qualidade;
-l <int>: tamanho da seed [5..*] (28);
-e <int>: total valor qualidade [10..30] máximo para as posições onde há mismatch, considerando o
alinhamento todo (70);
--maxbts <int> : número máximo de backtracks permitidos (125, 800 com --best);
--pairtries <int>: número máximo de tentativas de encontrar sequências em pares;
--try_hard: equivalente p/ valores altos de --maxbts e –pairtries
Reporting:
--best: reporta os melhores alinhamentos considerando número de mismatches na seed e o valor de
qualidade dessas bases; (1 mismatch qual 40 é melhor que 2 mismatches qual 10)
-a: reporta todos os alinhamentos válidos;
-k <int>: reporta até k alinhamentos válidos;
-m <int>: suprime os alinhamentos múltiplos de uma leitura se há mais que m alinhamentos válidos;
BWA
• http://bio-bwa.sourceforge.net/
• Dois algoritmos baseados na Transformação de
Burrows-Wheeler (BWT)
– sequências pequenas (queries) até ~200bp com baixa
taxa de erro (<3%)
• Alinhamento global (FM-index) com respeito às queries,
suporta sequências paired-ends.
• BWA (Li and Durbin, 2009)
– sequências longas (queries) com taxa de erro maior
• Alinhamento local = heurísticas+Smith-Waterman-like em
árvores de sufixos, não suporta sequências paired-ends.
Lento para sequências pequenas.
• BWA-SW (Li and Durbin, 2010)
BWA (index)
bwa index [-p prefix] [-a algoType] [-c] <in.db.fasta> Index database sequences in the
FASTA format.
Constrói os índices das sequências de referência que serão utilizadas no alinhamento.
OPTIONS:
-c
-p STR
-a STR
is
bwtsw
Build color-space index. The input fast should be in nucleotide
space.
Prefix of the output database [same as db filename]
Algorithm for constructing BWT index. Available options are:
IS linear-time algorithm for constructing suffix array.
- requires 5.37*N (N = size of the database)
- IS is moderately fast, but does not work with database larger
than 2GB.
Algorithm implemented in BWT-SW.
- does not work with database smaller than 10MB and it is usually
slower than IS.
[/data/indexes]$ bwa index –a bwtsw /data/hg19.fa
BWA (aln)
bwa aln [-n maxDiff] [-o maxGapO] [-e maxGapE] [-d nDelTail] [-i nIndelEnd] [-k maxSeedDiff] [-l seedLen] [-t nThrds] [-cRN] [-M misMsc] [-O gapOsc] [-E gapEsc] [-q
trimQual] [-f <out.sai>] <in.db.fasta> <in.query.fq>
Encontra as coordenadas das leituras no array de sufixos;
MAIN OPTIONS:
•
-I
The input is in the Illumina 1.3+ read format;
•
-n NUM
Maximum edit distance (maxDiff). Máximo de diferenças permitidas no alinhamento todo. Se inteiro é a distância de edição,
se ponto flutuante, automaticamente seleciona essa distância de edição com base nos valores de taxa de erro (uniforme) 2%,
tamanho das sequências e NUM (threshold) [0.04];
•
-l INT
Take the first INT subsequence as seed (seedLen) – Tamanho da primeira semente de alinhamento. Se INT maior que o tamanho
das sequências, desabilita seeding [32];
•
-k INT
Maximum edit distance in the seed (maxSeedDiff) – Número máximo de diferenças permitidas na primeira semente de
alinhamento (seed) [2];
•
-t INT
Number of threads (multi-threading mode) [1];
•
-o INT
Maximum number of gap opens [1]
•
-e INT
Maximum number of gap extensions, -1 for k-difference mode (disallowing long gaps) [-1]
•
-d INT
Disallow a long deletion within INT bp towards the 3’-end [16]
•
-i INT
Disallow an indel within INT bp towards the ends [5]
•
-N
Disable iterative search. All hits with no more than maxDiff differences will be found. This mode is much slower than the
default.
•
-q INT
Parameter for read trimming (soft clipping).
[/data/input]$ bwa aln –f out.sai /data/hg19.fa in.fastq
BWA (samse/sampe)
bwa samse [-n maxOcc] –f <out.sam> <in.db.fasta> <in.sai> <in.fq>
Gera alinhamentos no formato SAM a partir de leituras single-end. Hits repetitivos são escolhidos aleatoriamente.
[/data/input]$ bwa samse –f out.sam /data/hg19.fa out.sai in.fastq
bwa sampe [-a maxInsSize] [-o maxOcc] [-n maxHitPaired] [-N maxHitDis] [-P] -f <out.sam> <in.db.fasta> <in1.sai>
<in2.sai> <in1.fq> <in2.fq>
Gera alinhamentos no formato SAM a partir de leituras paired-end. Hits de pares repetitivos são escolhidos
aleatoriamente.
MAIN OPTIONS:
-n INT
Maximum number of alignments to output in the XA tag for reads paired properly. If a read
than INT hits, the XA tag will not be written. [3]
has more
-N INT
Maximum number of alignments to output in the XA tag for disconcordant read pairs
singletons). If a read has more than INT hits, the XA tag will not be written. [10]
(excluding
-s
Disable Smith-Waterman for the unmapped mate;
[/data/input]$ bwa sampe –f outpe.sam /data/hg19.fa out1.sai out2.sai \
> in1.fastq in2.fastq
BWA (bwasw)
bwa bwasw [-a matchScore] [-b mmPen] [-q gapOpenPen] [-r gapExtPen] [-t
nThreads] [-w bandWidth] [-T thres] [-s hspIntv] [-z zBest] [-N nHspRev] [-c
thresCoef] –f <out.sam> <in.db.fasta> <in.fq>
Gera alinhamento de sequências longas no formato SAM.
MAIN OPTIONS:
-t INT
Number of threads in the multi-threading mode [1]
-z INT
Z-best heuristics. Higher -z increases accuracy at the cost of
[1]
-s INT
Maximum SA interval size for initiating a seed. Higher -s
increases accuracy at the cost of speed. [3]
[/data/input]$ bwa bwasw –f out.sam /data/hg19.fa in.fastq
speed.
The Sequence Alignment Map
• Formato genérico para armazenar o resultado
dos alinhamentos de sequências de leituras
contra sequências de referência;
– (Li et al., 2009)
• Suporta o armazenamento de dados das mais
variadas plataformas de sequenciamento;
– Há duas seções:
• cabeçalho (iniciadas por @)
• corpo (alinhamentos);
– O alinhamento é representado no corpo por linhas com campos
(11 mandatórios) delimitados por TAB;
Formato SAM (1)
Col
1
2
3
4
5
6
7
8
9
10
11
12
Field
QNAME
FLAG
RNAME
POS
MAPQ
CIGAR
MRNM
MPOS
ISIZE
SEQ
QUAL
OPT
Description
Query (pair) NAME
bitwise FLAG
Reference sequence NAME
1-based leftmost POSition/coordinate of clipped sequence
MAPping Quality (Phred-scaled)
extended CIGAR string
Mate Reference sequence NaMe (‘=’ if same as RNAME)
1-based Mate POSistion
Inferred insert SIZE
query SEQuence on the same strand as the reference
query QUALity (ASCII-33 gives the Phred base quality)
variable OPTional fields in the format TAG:VTYPE:VALUE
Mapping Quality
• Qualidade do mapeamento – escala Phred de
qualidade (Li et al., 2008);
Formato SAM (2)
C
SAM
pile-up
forward match (.) mismatch (base in lower case)
reverse match (,) mismatch (base in upper case)
insertion +<LENGTH><BASES>
deletion -<LENGTH><BASES>
missing *
start read alignent ^
end read alignment $
Bitwise FLAG
Exemplo: FLAG 2 alinhamentos de leituras paired-end
FLAGs: 99 (00001100011 ) 147 (00010010011 )
Flag (Hex)
0x0001
0x0002
0x0004
0x0008
0x0010
0x0020
0x0040
0x0080
0x0100
0x0200
0x0400
Description
the read is paired in sequencing
the read is mapped in a proper pair
the query sequence itself is unmapped
the mate is unmapped
strand of the query (1 for reverse)
strand of the mate
the read is the first read in a pair
the read is the second read in a pair
the alignment is not primary
QC failure
optical or PCR duplicate
99
1
1
0
0
0
1
1
0
0
0
0
147
1
1
0
0
1
0
0
1
0
0
0
Formato BED
Formato flexível para definição de dados de anotação. Os três primerio campos são requeridos:
•
•
•
chrom
chromStart
chromEnd
– Nome do cromossomo
– Início da coordenada no cromossomo
– Final da coordenada no cromossomo
•
•
•
•
•
•
•
•
•
name
score
strand
thickStart
thickEnd
itemRgb
blockCount
blockSizes
blockStarts
– Nome da linha de anotação
– Pontuação entre 0 e 1000
– Fita do DNA + ou – Início da coordenada de destaque (por ex. start codon)
– Fim da coordenada de destaque (por ex. stop codon)
– Valor RGB no formato R,G,B (e.g. 255,0,0) ; cor no Genome Browser
– Número de blocos de alinhamento (por ex. exons)
– Lista (separada por vírgula) de tamanho dos blocos
– Lista (separada por vírgula) com o iníco da coordenada dos blocos
Exemplo
chr1 34610 36081 NR_026818 0 chr1 803450 812182 NR_027055 0 -
36081 36081 0 3 564,205,361, 0,666,1110,
812182 812182 0 3 605,1044,57, 0,6041,8675,
BEDtools
• http://code.google.com/p/bedtools/
• Conjunto de utilitários que permitem executar
tarefas comuns em genômica, tais como
encontrar coordenadas de alinhamentos que
possuem sobreposição;
• Funcionalidades
– bamToBed
Converte alinhamentos no formato
BAM para o formato BED
– intersectBed Retorna as sobreposições entre dois
arquivos no formato BED
– ...
samtools
• http://samtools.sourceforge.net/
• Provê várias ferramentas para manipulação dos alinhamentos no
formato SAM/BAM
• Funcionalidades
– import
– view
–
–
–
–
–
–
sort
merge
index
faidx
tview
pileup
Conversão de SAM-para-BAM
Conversão de BAM-para-SAM e recuperação de
subalinhamentos
Ordenação (posição nos cromossomos)
união de múltiplos alinhamentos ordenados
Indexação de alinhamentos ordenados
Indexação de arquivos FASTA e recuperação de subsequências
Visualização do alinhamento
Gera arquivo de cobertura por posição
goby
• http://campagnelab.org/software/goby/
• Sistema computacional para manipulação de dados da
nova geração de sequenciamento, com aplicações que
facilitam a implementação de “pipelines” de análise de
dados;
• Funcionalidades
– Filtragem de leituras redunantes;
– Criação de arquivos nos formatos wiggle/bed para
visualização;
– Análise de expressão gênica diferencial;
– Detecção de variações;
–…
Integrative Genomics Viewer
• IGV (Genome Browser)
– http://www.broadinstitute.org/software/igv/home
– Download
– ~10Mb (.zip)
– Permissão execução p/ todos
– chmod a+x igv_linux.sh
– Aumentar espaço de alocação de memória
– Editar arquivo igv_linux.sh (-Xmx1200m)
– Necessita acesso à internet!
Referências
•
•
•
•
•
•
•
•
Cock PJ, Fields CJ, Goto N, Heuer ML, Rice PM. The Sanger FASTQ file format for sequences with quality
scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res. 2010 Apr;38(6):1767-71. Epub 2009
Dec 16. Review. PubMed PMID: 20015970; PubMed Central PMCID: PMC2847217;
Li H, Homer N. A survey of sequence alignment algorithms for next-generation sequencing. Brief
Bioinform. 2010 Sep;11(5):473-83. Epub 2010 May 11. Review. PubMed PMID: 20460430; PubMed
Central PMCID: PMC2943993;
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform.
Bioinformatics. 2009 Jul 15;25(14):1754-60. Epub 2009 May 18. PubMed PMID: 19451168; PubMed
Central PMCID: PMC2705234;
Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics.
2010 Mar 1;26(5):589-95. Epub 2010 Jan 15. PubMed PMID: 20080505; PubMed Central PMCID:
PMC2828108;
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA
sequences to the human genome. Genome Biol. 2009;10(3):R25. Epub 2009 Mar 4. PubMed PMID:
19261174; PubMed Central PMCID: PMC2690996;
Li H, Ruan J, Durbin R. Mapping short DNA sequencing reads and calling variants using mapping quality
scores. Genome Res. 2008 Nov;18(11):1851-8. Epub 2008 Aug 19. PubMed PMID: 18714091; PubMed
Central PMCID: PMC2577856;
Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP. Integrative
genomics viewer. Nat Biotechnol. 2011 Jan;29(1):24-6. PubMed PMID: 21221095;
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome
Project Data Processing Subgroup. The Sequence Alignment/Map format and SAMtools. Bioinformatics.
2009 Aug 15;25(16):2078-9. Epub 2009 Jun 8. PubMed PMID: 19505943; PubMed Central PMCID:
PMC2723002;
Daniel Guariz Pinheiro
[email protected]