aula2_ibm1029

Prévia do material em texto

1
IBM1029
Introdução à Bioinformática
Aula 2
Profª Drª Silvana Giuliatti
Departamento de Genética – Bloco G
Ramal: 4503
silvana@rge.fmrp.usp.br
Faculdade de Medicina
Departamento de Genética
PHRED
Genome Research 8: 175-185, 1998
2
Genome Research 8: 186-194, 1998
¾Lê dados de sinais de fluorescência vindos de 
sequenciadores automáticos de DNA.
¾ Define as bases (base calling)
¾ Atribui valores de qualidade às bases: um valor 
baseado na estimação do erro calculado para cada base 
individualmente. 
¾ Gera arquivo de saída: bases definidas e valores de 
qualidade.
Traço Ideal: picos sem sobreposição, igualmente espaçados.
Traço real: os picos desviam da posição ideal por 
imperfeições no sequenciamento. 
O algoritmo pode ser dividido em fases:
1) Determinação das posições dos picos ideais. 
2) Determinação das posições dos picos observados. 
3)Comparação entre os picos ideais e observados através de 
programação dinâmica.
4) Quantificação da qualidade do sinal.
3
Fase 1: Determinação das posições dos picos ideais.
• Baseado no fato de que os fragmentos são, 
relativamente, igualmente espaçados, em média, na 
maioria das regiões do gel.
• A predição dos picos ideais é feita usando Fourier. 
Inicia o processo de predição a partir da análise 
dos quatro traços.
Fase 2: Determinação das posições dos picos 
observados.
• Phred examina cada traço para encontrar o centro dos 
picos reais e suas áreas relativas aos seus vizinhos.
• Os picos são detectados independentemente ao longo de 
cada traço, uma vez que pode haver picos sobrepostos.
Fase 3: Comparação entre os picos ideais e 
observados através de programação dinâmica.
• O algoritmo de programação dinâmica é usado para 
comparar os picos reais detectados na fase 2 com os picos 
ideais localizados na fase 1.
Fase 4: Quantificação da qualidade do sinal.
a) Suponha:
B = 1.000.000 bases
E = 10.000 erros
P = 10.000/1.000.000 = 0,01 (probabilidade de erro)
b) Suponha: B/2 → B1 = B2 = 50.000 (sub-conjuntos de bases) e 
E/2 → E1 = 9.000 erros e E2 = 1.000 erros
P1 = 9.000/500.000 = 0,018
P2 = 1.000/500.000 = 0,002
P1 + P2 = 0,02 (probabilidade de erro)
Melhor maneira de discriminar as bases mais exatas das menos exatas.
4
Objetivo: desenvolver probabilidades de erro que sejam válidas 
e que tenham alto poder de discriminação para baixos 
valores de erro (r).
• Algoritmo novo para derivar as probabilidades de erro.
Parâmetros Utilizados
1. Espaçamento de Picos: janela com 7 picos.
2. Razão Base não-identificada/Base identificada:
janela com 7 picos.
3. Mesmo que item 2, mas com janela de 3 picos.
4. Resolução do Pico
• Valor limite para cada parâmetro: ri, sj, tk, um
Definições
err(i, j, k, m) : número total de bases erradas abaixo do corte que 
satisfaça: r(b) ≤ ri, s(b) ≤ sj, t(b) ≤ tk, u(b) ≤ um
corr(i, j, k, m) : número total de bases corretas abaixo do corte.
Taxa de Erro abaixo do corte: p = (1+ err)/(1+err+corr)
Onde, o valor 1 é um valor de correção adicionado para ter 
certeza de que ambos, numerador e denominador, serão 
positivos.
Valor de Qualidade
q = - 10 x log (p)
onde:
p - probabilidade do erro estimado para aquela base-call
Exemplo:
q = 20 significa p = 10-2 (1 erro em 100 bases)
Valores de qualidades mais altos, representam 
probabilidades de erro menores
5
Execução 
• Opções de Entrada
• -id <nome_dir>: lê e processa arquivo no diretório 
nome_dir.
• -if <nome_arq>: lê e processa arquivos listados no 
arquivo nome_arq.
• Opções de Processamento
• -trim-alt <seq>: faz um corte na seqüência. Localiza 
a região de alta qualidade de cada read. Para cada base, 
a probabilidade de erro é subtraído do valor default ou 
do valor selecionado em -trim-cutoff.
• -trim-cutoff <valor>: seleciona o valor para a 
probabilidade de erro. O default é 0.05.
• Opções de Saída
• -pd <nome_dir>: escreve arquivos PHD no diretório 
<nome_dir>.
• -sa <nome_arq>: escreve uma sequência de saída no 
formato FASTA com o nome nome_arq.
• -qa <nome_arq>: escreve arquivo de saída das 
qualidades num formato FASTA com o nome 
nome_arq.
[root@localhost edit_dir]# phred -id ../chromat_dir -pd
../phd_dir
../chromat_dir/CM0-TF01020-C28F
../chromat_dir/CM0-TF01020-C28R
[root@localhost aula2]# cd phd_dir/
[root@localhost phd_dir]# ls
CM0-TF01020-C28F.phd.1 CM0-TF01020-C28R.phd.1
6
A estrutura de Arquivo phd
BEGIN_SEQUENCE CM0-TF01020-C28F
BEGIN_COMMENT
CHROMAT_FILE: CM0-TF01020-C28F
ABI_THUMBPRINT: 
067262134022000116163260066255
PHRED_VERSION: 0.990722.f
CALL_METHOD: phred
QUALITY_LEVELS: 99
TIME: Sun Feb 20 17:47:29 2005
TRACE_ARRAY_MIN_INDEX: 0
TRACE_ARRAY_MAX_INDEX: 10951
TRIM: 18 274 0.0500
CHEM: term
DYE: big
END_COMMENT
BEGIN_DNA
c 4 6
c 4 23
g 6 36
t 4 38
a 4 55
t 4 68
g 6 79
g 4 93
c 4 112
g 4 117
g 4 134
t 4 137
t 4 160
n 0 169
g 4 171
a 9 192
a 16 206
a 9 220
a 18 233
c 18 243
g 32 253
a 29 266
c 37 278
...
c 11 10755
a 9 10766
c 4 10781
g 4 10790
t 9 10798
a 9 10808
t 10 10821
g 6 10831
t 6 10840
g 6 10849
c 6 10850
t 8 10865
c 8 10882
t 8 10890
c 4 10908
a 4 10911
c 4 10927
c 6 10937
c 6 10939
END_DNA
END_SEQUENCE
g 4 10522
a 4 10529
g 4 10543
g 4 10551
t 8 10567
t 8 10579
a 6 10589
c 6 10593
g 11 10603
t 6 10622
t 8 10631
a 4 10643
g 4 10653
c 8 10668
g 6 10678
g 6 10682
g 8 10698
c 6 10706
a 6 10707
a 8 10722
c 9 10737
Cabeçalho
Base
Valor de 
qualidade
Posição da base 
no traço
Phd2Fasta
[root@localhost edit_dir]# phd2fasta -id ../phd_dir -os fasta_seq -oq
fasta_seq.screen.qual 
[root@localhost edit_dir]# ls
fasta_seq fasta_seq.screen.qual
[root@localhost edit_dir]#more fasta_seq
>CM0-TF01020-C28F CHROMAT_FILE: CM0-TF01020-C28F PHD_FILE: CM0-
TF01020-C28F.phd.1 CHEM: term DYE: big TIME: Sun Mar 6 00:37:30 2005
ccgtatggcggttngaaaacgacggccagtgccaagcttgcatgcctgca
ggtcgactctagaggatcccctgactgggtttgttccatcccagtcagcg
atgtttatgtggccctgactgggtttgttccatccgatgatgattccatt
cgattccgttcaatgattattccattcgagtccattcgatgattccattt
gattccattcgatgatgattgcattcgagtccatggattattccattcca
ttccattagatgattccattcgggtccgttcgaagattctcttcgattcc
attcgatgattccgtttttttccgtttggtgttgataccattcgatttcc
attcgatgataattccattcgatctatgcgagntttaattggggcccggg
ccacataaacatcgccatcatcaggaaatgcactcanaacacagtgagat
ccatctacccagtagagtgcantattaaaagacagaatacagattctgtg
agagcagnaagatatgatatcttttgtggaatgtaattngatactctatg
more fasta_seq.screen.qual
>CM0-TF01020-C28F PHD_FILE: CM0-TF01020-C28F.phd.1
4 4 6 4 4 4 6 4 4 4 4 4 4 0 4 9 16 9 18 18 32 29 37
34 26 26 23 26 22 25 22 20 20 22 22 25 29 39 37 37
37 33 31 34 32 29 29 27 27 29 29 32 34 40 33 32 34
36 39 39 29 29 29 29 29 35 40 40 40 40 56 40 40 37
35 35 34 32 34 34 35 35 37 37 37 39 40 40 45 51 46
46 42 37 44 44 48 56 56 56 56 56 51 42 46 42 42 42
44 35 35 35 35 35 35 35 35 35 32 29 29 29 37 37 37
35 28 28 28 35 35 35 37 37 32 29 26 26 26 26 26 24
33 31 35 36 35 37 42 42 36 35 32 32 32 23 16 11 11
19 24 31 35 35 35 30 30 35 35 35 35 35 29 29 35 35
35 33 32 23 23 23 29 32 35 37 37 37 40 42 42 35 35
35 35 35 35 40 45 40 40 37 37 35 35 32 35 38 40 40
40 44 44 44 44 44 44 44 33 33 21 18 18 33 33 
7
Cross_Match
Retira as bases contaminantes que possam atrapalhar na 
montagem pelo Phrap. 
Características:
88 Compara cada seqüência no formato FASTA, gerada 
pelo programa phd2fasta, com as seqüências contidas no 
arquivo vector.seq. 
8 Baseado numa implementação do algoritmo de 
Smith-Waterman-Gotoh. 
8 “Mascara” as bases contaminantes substituindo-as por 
“X”.
88 A seqüência com máscara é colocada num arquivo 
com extensão .screen.
8 É executado automaticamente quando o script 
phredphrap é executado. 
[root@localhost edit_dir]# cross_match fasta_seq -minmatch 20 -screen
> screen.out
cross_match fasta_seq -minmatch 20 -screencross_match version 0.990329
Reading parameters ... 0.504 Mbytes allocated -- total 0.504 Mbytes
Run date:time 050308:085140
Done
Total space allocated: 0.504 Mbytes; currently free: 0.499 Mbytes in 
2 blocks
Reading query file into memory ... Done
Total space allocated: 0.504 Mbytes; currently free: 0.181 Mbytes in 
2 blocks
Complementing ... Done
Total space allocated: 0.504 Mbytes; currently free: 0.181 Mbytes in 
2 blocks
Allocating align_entries ... Done
Total space allocated: 0.504 Mbytes; currently free: 0.055 Mbytes in 
2 blocks
Reading quality files ...
...
8
[root@localhost edit_dir]# ls
fasta_seq fasta_seq.log fasta_seq.screen fasta_seq.screen.qual 
screen.out
[root@localhost edit_dir]# more fasta_seq.screen
>CM0-TF01020-C02F CHROMAT_FILE: CM0-TF01020-C02F PHD_FILE: CM0-
TF01020-C02F.phd.1 CHEM: term DYE: big TIME: Sun Mar 6 10:20:47 
2005
TCXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXGCGGGAGTAAACTTGCCCGTGGAAGGAGTTTGGGCGGA
CCGCTATTCAACGGCAAGTTTTTTGGAAGAGNANNNNNNNNNNNNNANTT
NNNNNNNNNNNNNNNNNNNNTTGGTTGGTGGGGTGGTTGNNGNGGGTGTG
TNNGGTGTGTGTGTGTTGTTGGGTTTTTTTGTTGGNTGTNTTTNTGGTTT
TGTTTTTGTTGGTGTTTGTTGGTGTGTGGGTTGTGGTGGGGGTGTTGGTT ...
[root@localhost edit_dir]# more screen.out
cross_match fasta_seq -minmatch 20 -screen
cross_match version 0.990329
Run date:time 050308:085140
Query file(s): fasta_seq
Presumed sequence type: DNA
Pairwise comparison algorithm: banded Smith-Waterman
Score matrix (set by value of penalty: -2)
A C G T N X
A 1 -2 -2 -2 0 -3
C -2 1 -2 -2 0 -3
G -2 -2 1 -2 0 -3
T -2 -2 -2 1 0 -3
N 0 0 0 0 0 0
X -3 -3 -3 -3 0 -3
Gap penalties: gap_init: -4, gap_ext: -3, ins_gap_ext: -3, 
del_gap_ext: -3,
Using complexity-adjusted scores. Assumed background frequencies:
A: 0.250 C: 0.250 G: 0.250 T: 0.250 N: 0.000 X: 0.000
...