Biologia universitária — 3.º ano · Bachelor Year 3

5Bioinformática e análise de sequências

Um biólogo que acaba de sequenciar um gene de um verme do mar profundo cola seus $300$ aminoácidos num formulário da web e, três segundos depois, fica sabendo que a proteína é uma prima distante de uma cinase humana, com $31\,\%$ de identidade ao longo de $280$ resíduos e uma probabilidade de $10^{-40}$ de a semelhança ser acaso. Por trás desses três segundos estão um algoritmo de programação dinâmica de 1970, uma teoria estatística dos alinhamentos aleatórios, matrizes de substituição destiladas de milhares de famílias de proteínas e um banco de dados de umas cem bilhões de posições. Este capítulo trata do raciocínio dentro da caixa: como duas sequências são alinhadas de modo que o alinhamento seja comprovadamente o melhor, como o escore passa a significar alguma coisa, como se distingue uma correspondência de uma coincidência e como se acham padrões num genoma que ninguém examinou antes. A matemática é elementar — uma recorrência, um logaritmo, uma distribuição de Poisson — e vale a pena conhecê-la, porque toda conclusão tirada de uma comparação de sequências repousa nela.

5.1 Alinhar duas sequências

Definição 5.1 (Alinhamento e escore)

Um alinhamento de duas sequências as escreve uma sobre a outra, com lacunas (–) inseridas de modo que as colunas emparelhem um resíduo com um resíduo ou um resíduo com uma lacuna, e nenhuma coluna emparelhe duas lacunas. Seu escore é a soma, sobre as colunas, de um escore de substituição $s(a,b)$ para cada par de resíduos e de uma penalidade de lacuna para cada lacuna: uma penalidade linear $-d$ por posição de lacuna ou, mais realisticamente, uma penalidade afim $-d - (k-1)e$ para uma série de $k$ lacunas, com o custo de abertura $d$ maior que o custo de extensão $e$ , já que uma inserção de vários resíduos é um único evento evolutivo. Um alinhamento global cobre as duas sequências de ponta a ponta; um alinhamento local encontra o par de subcadeias de maior escore e ignora o resto, que é o que se quer quando um domínio compartilhado está em duas proteínas de resto não aparentadas.

Teorema 5.2 (Needleman–Wunsch)

Sejam $x = x_{1}\dots x_{m}$ e $y = y_{1}\dots y_{n}$ , com penalidade linear de lacuna $d$ . Defina $F(i,j)$ como o melhor escore de um alinhamento global dos prefixos $x_{1}\dots x_{i}$ e $y_{1}\dots y_{j}$ . Então $F(i,0) = -id$ , $F(0,j) = -jd$ e, para $i,j \ge 1$ ,

F(i,j) = \max\bigl\{\,F(i-1,j-1) + s(x_{i},y_{j}),\; F(i-1,j) - d,\; F(i,j-1) - d\,\bigr\}.

$F(m,n)$ é o escore global ótimo, um alinhamento ótimo se recupera retrotraçando a partir de $(m,n)$ as escolhas que produziram cada máximo, e todo o cálculo leva $mn$ passos. A variante Smith–Waterman para alinhamento local acrescenta $0$ como quarta opção no máximo, fixa as bordas em $0$ e lê a resposta na maior entrada da tabela.

Demonstração. Considere a última coluna de qualquer alinhamento dos dois prefixos. Ela é uma de três coisas: $x_{i}$ sobre $y_{j}$ , $x_{i}$ sobre uma lacuna, ou uma lacuna sobre $y_{j}$ . Removê-la deixa um alinhamento de $(x_{1}\dots x_{i-1}, y_{1}\dots y_{j-1})$ , de $(x_{1}\dots x_{i-1}, y_{1}\dots y_{j})$ ou de $(x_{1}\dots x_{i}, y_{1}\dots y_{j-1})$ , respectivamente, cujo escore é no máximo o $F$ desse par; e, ao contrário, cada um desses alinhamentos ótimos pode ser estendido pela coluna final correspondente. Logo, o melhor escore que termina em cada tipo de coluna é o $F$ do par mais curto mais o escore da coluna, e o ótimo é o maior dos três. As bordas são forçadas (só lacunas são possíveis contra um prefixo vazio). A indução sobre $i + j$ preenche a tabela; o número de células é $(m+1)(n+1)$ . Para o alinhamento local, a opção extra $0$ significa “começar um novo alinhamento aqui”, o que faz de $F(i,j)$ o melhor escore de um alinhamento que termina em $(i,j)$ , e o melhor alinhamento local termina em algum lugar. ∎

Exemplo 5.3 (Uma tabela quatro por três)

Alinhe GAT com GCAT, pontuando $+1$ para correspondência, $-1$ para discordância, $d = 1$ . As bordas são $0, -1, -2, -3, -4$ ao longo do topo e $0, -1, -2, -3$ pela lateral. Preenchendo linha a linha: $F(\text{G},\text{G}) = 1$ , $F(\text{G},\text{C}) = 0$ , $F(\text{G},\text{A}) = -1$ , $F(\text{G},\text{T}) = -2$ ; $F(\text{A},\text{G}) = 0$ , $F(\text{A},\text{C}) = 0$ , $F(\text{A},\text{A}) = 1$ , $F(\text{A},\text{T}) = 0$ ; $F(\text{T},\text{G}) = -1$ , $F(\text{T},\text{C}) = -1$ , $F(\text{T},\text{A}) = 0$ , $F(\text{T},\text{T}) = 2$ . O ótimo é $2$ e o retrotraçado — diagonal a partir de (T,T), diagonal a partir de (A,A), depois à esquerda de (G,C) até (G,G), depois diagonal — dá

\begin{array}{c} \texttt{G-AT}\\ \texttt{GCAT} \end{array}

três correspondências e uma lacuna: $3 - 1 = 2$ .

A tabela de Needleman–Wunsch para GAT contra GCAT (correspondência +1, discordância -1, lacuna -1). Cada célula é o melhor escore para os dois prefixos que ali terminam; o caminho vermelho, retrotraçado a partir do canto, é o alinhamento ótimo. — A tabela de Needleman–Wunsch para GAT contra GCAT (correspondência $+1$ , discordância $-1$ , lacuna $-1$ ). Cada célula é o melhor escore para os dois prefixos que ali terminam; o caminho vermelho, retrotraçado a partir do canto, é o alinhamento ótimo.

Método 5.4 (Alinhar duas sequências)

(1) Escolha a pontuação: uma matriz de substituição adequada à divergência esperada (BLOSUM62 para proteínas de distância desconhecida; correspondência/discordância para DNA) e penalidades afins de lacuna (tipicamente abertura $-11$ , extensão $-1$ com a BLOSUM62). (2) Decida entre global e local: global para duas sequências que se acreditam homólogas em todo o seu comprimento, local nos demais casos. (3) Preencha a tabela pela recorrência, guardando para cada célula um ponteiro para a escolha que deu seu máximo. (4) Retrotrace a partir de $(m,n)$ (global) ou da célula máxima até um zero (local), escrevendo o alinhamento da direita para a esquerda. (5) Julgue o resultado não por seu escore bruto, mas por sua significância estatística (adiante), e olhe para ele: lacunas longas, trechos de baixa complexidade e alinhamentos confinados a uma repetição são avisos.

5.2 Pontuar: quanto vale uma correspondência

Definição 5.5 (Matrizes de substituição)

Uma matriz de substituição dá $s(a,b)$ para cada par de aminoácidos como um escore de log-chances:

s(a,b) = \frac{1}{\lambda}\,\log\frac{q_{ab}}{p_{a}\,p_{b}},

em que $q_{ab}$ é a frequência com que $a$ e $b$ são encontrados alinhados em alinhamentos confiáveis de proteínas aparentadas, $p_{a} p_{b}$ é a frequência com que seriam emparelhados por acaso, e $\lambda$ é uma escala escolhida para que as entradas sejam inteiros convenientes. Um escore positivo significa que o par ocorre mais vezes em homólogas do que por acaso; os escores de identidade são maiores para os aminoácidos raros (triptofano $+11$ , cisteína $+9$ na BLOSUM62) e menores para os comuns (leucina $+4$ , alanina $+4$ ), e as substituições conservativas (isoleucina–valina $+3$ ) pontuam positivo enquanto as radicais (triptofano–glicina $-2$ ) pontuam negativo. As matrizes PAM (Dayhoff, 1978) foram derivadas de proteínas proximamente aparentadas e extrapoladas para distâncias maiores por multiplicação de matrizes; as matrizes BLOSUM (Henikoff e Henikoff, 1992) foram contadas diretamente em blocos de sequências alinhadas agrupadas a uma dada identidade — a BLOSUM62 a partir de blocos a $62\,\%$ — e são o padrão porque foram medidas, e não extrapoladas, na distância em que são usadas.

Proposição 5.6 (Por que log-chances)

Para que um esquema de pontuação possa ser usado em alinhamento local, o escore esperado de uma coluna emparelhada ao acaso, $\sum_{a,b} p_{a} p_{b}\, s(a,b)$ , precisa ser negativo, e alguns escores precisam ser positivos; do contrário os alinhamentos aleatórios cresceriam sem limite e o segmento de maior escore seria a sequência inteira. Dado isso, qualquer esquema assim é equivalente a um esquema de log-chances para algumas frequências-alvo $q_{ab}$ — os alinhamentos que ele achará ótimos são aqueles cujos pares de resíduos se distribuem como $q_{ab}$ . Escolher a matriz é, portanto, escolher a divergência que se espera detectar: uma matriz para parentes próximos (BLOSUM80, PAM30) tem positivos mais agudos e negativos mais duros, e uma para parentes distantes (BLOSUM45, PAM250) é mais achatada.

Demonstração. Admitido neste nível. ∎

Exemplo 5.7 (Identidade, similaridade e a zona crepuscular)

Duas sequências proteicas aleatórias alinhadas de forma ótima com lacunas chegam a cerca de $15\text{ a }20\,\%$ de identidade por acaso. Acima de $35\,\%$ de identidade ao longo de cem resíduos, duas proteínas são quase certamente homólogas; entre $20\,\%$ e $35\,\%$ está a zona crepuscular, em que a identidade sozinha não decide e a estatística adiante precisa decidir. As homólogas podem cair bem abaixo da zona: as subunidades da hemoglobina e a mioglobina compartilham $25\,\%$ de identidade, a lisozima e a $\alpha$ -lactalbumina $40\,\%$ , e muitos pares de proteínas de mesmo enovelamento compartilham menos de $15\,\%$ , detectáveis só pela comparação de perfis ou de estruturas.

5.3 Buscar num banco de dados

Definição 5.8 (BLAST)

Alinhar uma consulta de $300$ resíduos contra um banco de $10^{11}$ por programação dinâmica completa custaria $3\times 10^{13}$ atualizações de célula por busca. O BLAST (Altschul e colaboradores, 1990) troca um pouco de sensibilidade por mil vezes mais velocidade em três etapas: (1) lista as palavras da consulta (três resíduos para proteínas, onze bases para DNA) e suas vizinhas de alto escore; (2) varre o banco em busca de correspondências exatas de palavra — as sementes; (3) estende cada semente nos dois sentidos sem lacunas até que o escore caia uma quantidade fixada abaixo de seu melhor valor, guardando os pares de segmentos de alto escore (HSPs) e unindo depois os HSPs próximos por programação dinâmica com lacunas numa faixa estreita. Uma homóloga verdadeira quase sempre contém ao menos uma palavra exata de três resíduos em comum; uma semelhança de acaso raramente contém, e nunca é estendida.

A heurística do BLAST. Palavras exatas e curtas compartilhadas pela consulta e pela entrada do banco (vermelho) são sementes; cada uma é estendida ao longo de sua diagonal enquanto o escore continua subindo, e só as extensões que se mantêm altas se tornam pares de segmentos de alto escore.

Teorema 5.9 (A estatística de um acerto ao acaso)

Para uma consulta de comprimento $m$ buscada contra um banco de comprimento total $n$ , com um esquema de pontuação de escore esperado negativo, o número de alinhamentos locais sem lacunas que pontuam ao menos $S$ e surgem por acaso segue uma distribuição de Poisson de média

E = K\,m\,n\,e^{-\lambda S},

em que $\lambda$ e $K$ dependem apenas do esquema de pontuação e das frequências dos resíduos ( $\lambda$ é a escala da matriz de log-chances). $E$ é o valor esperado do escore $S$ . Escrevendo o escore em bits, $S' = (\lambda S - \ln K)/\ln 2$ , a fórmula se torna $E = m n\, 2^{-S'}$ , e a probabilidade de ao menos um alinhamento de acaso alcançar $S$ é $P = 1 - e^{-E}$ , que é igual a $E$ quando $E$ é pequeno.

Prova parcial. A cauda exponencial é o teorema de Karlin–Altschul e é admitida: o escore máximo de segmento de um passeio aleatório com deriva negativa tem uma distribuição cuja cauda decai como $e^{-\lambda S}$ , com $\lambda$ a raiz positiva de $\sum_{a,b} p_{a} p_{b} e^{\lambda s(a,b)} = 1$ — que é exatamente a equação que torna a matriz de log-chances consistente. Dada essa cauda, o resto é contagem. Os segmentos de alto escore podem começar em qualquer um dos $mn$ pares de posições, são raros e são quase independentes; o número deles que excede $S$ é, portanto, Poisson com média proporcional a $mn$ e à probabilidade da cauda, $E = Kmn\,e^{-\lambda S}$ . A probabilidade de nenhum é $e^{-E}$ . A substituição pelo escore em bits é álgebra: $e^{-\lambda S} K = 2^{-(\lambda S - \ln K)/\ln 2}$ . Para alinhamentos com lacunas vale a mesma forma, com $\lambda$ e $K$ estimados por simulação. ∎

Exemplo 5.10 (Ler um valor E)

Uma consulta de $250$ resíduos contra um banco de $5\times 10^{10}$ resíduos tem $mn = 1.25\times 10^{13} \approx 2^{43.5}$ . Um acerto com escore em bits de $60$ tem $E = 2^{43.5 - 60} = 2^{-16.5} \approx 10^{-5}$ : é, essencialmente com certeza, uma homóloga. Um acerto com $S' = 40$ tem $E = 2^{3.5} \approx 11$ : esperam-se onze escores desses por acaso, e o acerto não significa nada. O mesmo alinhamento, com o mesmo escore em bits, buscado contra um banco dez vezes maior, tem um $E$ dez vezes maior — a significância é propriedade da busca, não do par. O limiar de uso comum é $E < 10^{-3}$ para uma homóloga confiável; $E \approx 0.01$ – $1$ merece uma segunda olhada com um método de perfil.

E = mn\,2-S': cada bit a mais reduz à metade o número esperado de acertos ao acaso, e um banco dez vezes maior custa 3.3 bits de significância para o mesmo alinhamento. — $E = mn\,2^{-S'}$ : cada bit a mais reduz à metade o número esperado de acertos ao acaso, e um banco dez vezes maior custa $3.3$ bits de significância para o mesmo alinhamento.

5.4 Perfis, estados ocultos e motivos

Definição 5.11 (Alinhamento múltiplo e perfis)

Um alinhamento múltiplo de sequências organiza uma família de sequências em colunas de resíduos homólogos. A programação dinâmica exata sobre $k$ sequências custa $n^{k}$ e é impossível além de três; os programas práticos alinham progressivamente, primeiro o par mais próximo segundo uma árvore-guia, depois sequências e grupos ao alinhamento em crescimento, com rodadas de refinamento. Um alinhamento pronto é resumido como um perfil: para cada coluna, a frequência de cada resíduo e das lacunas. Um modelo oculto de Markov de perfil formaliza isso como uma cadeia de estados de correspondência, um por coluna conservada, cada um emitindo resíduos com suas próprias probabilidades, com estados de inserção e de deleção que admitem resíduos a mais ou a menos em cada posição; o modelo de uma família (uma entrada do Pfam) pontua uma nova sequência pela probabilidade do melhor caminho pelos estados, e encontra homólogas bem abaixo da zona crepuscular da comparação par a par, porque uma coluna que só tolera resíduos hidrofóbicos diz isso, ao passo que uma sequência isolada não pode dizer.

Um modelo oculto de Markov de perfil de uma família de quatro colunas. Cada estado de correspondência M emite um resíduo com as frequências próprias da coluna; os estados de inserção I (com laços sobre si) admitem resíduos extras, e os estados de deleção D saltam uma coluna. Pontuar uma sequência é encontrar seu caminho mais provável.

Definição 5.12 (Motivos e conteúdo de informação)

Um motivo é um padrão curto — um sítio de fator de transcrição, um sinal de splicing, um sítio de fosforilação — representado por uma matriz de pesos por posição com a frequência $f_{i}(b)$ de cada base ou resíduo $b$ em cada posição $i$ . O conteúdo de informação da posição $i$ é $R_{i} = 2 - H_{i}$ bits para o DNA, em que $H_{i} = -\sum_{b} f_{i}(b)\log_{2} f_{i}(b)$ é sua entropia: $2$ bits para uma base invariante, $0$ para uma posição em que as quatro são igualmente prováveis. O total $R = \sum_{i} R_{i}$ é desenhado como um logo de sequência, cada posição uma pilha de letras cuja altura total é $R_{i}$ e cujas letras têm tamanho proporcional à frequência.

Proposição 5.13 (Quanta informação um sítio precisa ter)

Um sítio que precisa ser encontrado $\gamma$ vezes num genoma de $G$ posições, e em nenhum outro lugar, precisa de cerca de $R_{\text{necessário}} = \log_{2}(G/\gamma)$ bits de conteúdo de informação: o motivo precisa reduzir as $G$ posições candidatas às $\gamma$ verdadeiras, e cada bit reduz as candidatas à metade. Os motivos observados de reguladores bacterianos bem estudados correspondem a essa previsão — os sítios de E. coli de um repressor que se liga a algumas dezenas de lugares num genoma de $4.6\,\mathrm{Mb}$ carregam $16\text{ a }18$ bits; os motivos de fatores de transcrição eucariontes, com $8\text{ a }12$ bits num genoma de $3\times 10^{9}$ , não conseguem especificar sozinhos seus alvos, e é por isso que agem em combinações e na cromatina aberta do Capítulo 1.

Demonstração. Uma posição aleatória corresponde a um motivo de conteúdo de informação $R$ com probabilidade de cerca de $2^{-R}$ (cada bit de especificidade reduz a chance à metade), de modo que o número esperado de correspondências ao acaso em $G$ posições é $G\,2^{-R}$ . Para que os sítios verdadeiros se destaquem, isso precisa ser da ordem de $\gamma$ ou menos: $G\,2^{-R} \le \gamma$ , isto é, $R \ge \log_{2}(G/\gamma)$ . ∎

Exemplo 5.14 (Correspondências esperadas ao acaso)

Um sítio de restrição de seis bases fixas tem $R = 12$ bits e corresponde a uma posição aleatória com probabilidade $4^{-6} = 2^{-12}$ : cerca de $1100$ vezes num genoma de E. coli de $4.6\,\mathrm{Mb}$ lido nas duas fitas (o sítio é palindrômico, e portanto uma vez por posição) e $7\times 10^{5}$ vezes no genoma humano. Um fator eucarionte cujo motivo carrega $10$ bits corresponde a $3\times 10^{9}\times 2^{-10} \approx 3$ milhões de posições no genoma humano, vários milhares de vezes mais do que os genes que ele regula. Um motivo sozinho é um previsor fraco num genoma grande; o estado da cromatina, os motivos vizinhos e a conservação do sítio entre espécies é que fazem uma previsão.

Um logo de sequência de um motivo promotor do tipo caixa TATA. A altura de cada pilha é o conteúdo de informação daquela posição,

2 - H_{i}

bits; as quatro primeiras posições são quase invariantes e carregam a maior parte dos

12

bits do motivo.

5.5 Da sequência à função

Método 5.15 (Anotar uma proteína desconhecida)

Dada uma nova sequência codificadora: (1) traduza-a na fase certa e procure um peptídeo-sinal, segmentos transmembrana e regiões de baixa complexidade; (2) busque nos bancos de proteínas com o BLAST e leia os acertos com $E < 10^{-3}$ , observando se o alinhamento cobre a proteína inteira (uma ortóloga verdadeira) ou um segmento (um domínio compartilhado); (3) busque nos bancos de domínios com HMMs de perfil, que encontram famílias que o BLAST deixa passar e repartem a proteína em domínios; (4) infira ortologia, e não mera similaridade, verificando se o melhor acerto no outro genoma tem a consulta como seu melhor acerto (melhores acertos recíprocos) ou situando a proteína numa árvore de genes (Capítulo 25); (5) transfira a função das ortólogas com cautela — um resíduo catalítico conservado argumenta a favor de uma química conservada, e um resíduo ausente argumenta contra — e preveja a estrutura; (6) trate toda previsão como hipótese para a bancada.

Proposição 5.16 (Estrutura a partir da sequência)

O enovelamento de uma proteína é determinado por sua sequência (Capítulo 7), e calculá-lo a partir da sequência foi durante cinquenta anos o problema central não resolvido da área. Três abordagens tiveram êxito, uma após a outra. A modelagem por homologia constrói a estrutura de uma proteína sobre a de uma homóloga já resolvida, de modo confiável acima de $30\,\%$ de identidade. A análise de coevolução explora o fato de que dois resíduos em contato no enovelamento tendem a mutar juntos ao longo de um alinhamento múltiplo profundo, de modo que pares de colunas estatisticamente acopladas são contatos previstos, e contatos bastantes definem um enovelamento. Os métodos de aprendizado profundo treinados nas cem mil estruturas resolvidas e nesses alinhamentos agora preveem a maioria das estruturas de proteínas globulares com acurácia quase experimental (as avaliações CASP de 2020), e há bancos com uma estrutura prevista para praticamente toda sequência proteica conhecida. O que eles preveem pior é aquilo que uma estrutura única não capta: regiões desordenadas, conformações alternativas, o efeito de uma mutação pontual e os complexos.

À esquerda: uma estrutura proteica prevista, colorida pela confiança do modelo, de alta (azul) a baixa (laranja) numa alça desordenada. À direita: um escritório de bioinformática — navegadores de genomas e árvores nas telas, e nenhuma bancada úmida à vista.

Observação 5.17 (Os limites da inferência)

A maioria das anotações funcionais nos bancos nunca foi testada; elas foram transferidas de uma homóloga, que por sua vez fora anotada por transferência. Os erros se propagam e se multiplicam, e uma anotação errada numa proteína bem conectada pode infectar toda uma família. Os remédios são os de cima: distinga ortologia de homologia, leia o alinhamento, procure os resíduos catalíticos e lembre que “proteína hipotética” é um rótulo honesto que um terço dos genes da maioria dos genomas ainda merece.

5.6 Exercícios

Exercício 5.1 ★

Defina alinhamento global e local e dê uma situação biológica que exija cada um.

Solução

Solução de Exercício 5.1.

Global: as duas sequências alinhadas de ponta a ponta, cada resíduo numa coluna — para duas proteínas que se acreditam homólogas em todo o comprimento, como as ortólogas de uma enzima de manutenção. Local: o par de subcadeias de maior escore, ignorado o resto — para achar um domínio compartilhado (um domínio SH2 em duas proteínas de sinalização de resto não aparentadas) ou um gene numa longa sequência genômica.

Exercício 5.2 ★

Preencha a tabela de Needleman–Wunsch para AGC contra AAC com correspondência $+1$ , discordância $-1$ , lacuna $-1$ , e dê o alinhamento ótimo e o escore.

Solução

Solução de Exercício 5.2.

Bordas $0,-1,-2,-3$ nos dois sentidos. Linha A: $1, 0, -1$ . Linha G: $0, 0, -1$ . Linha C: $-1, -1, 1$ . Ótimo $F(3,3) = 1$ : AGC sobre AAC sem lacunas (correspondência, discordância, correspondência: $1 - 1 + 1 = 1$ ).

Exercício 5.3 ★

Na BLOSUM62, triptofano–triptofano pontua $+11$ e leucina–leucina $+4$ . Explique, a partir da fórmula de log-chances, por que a identidade do resíduo mais raro vale mais.

Solução

Solução de Exercício 5.3.

$s(a,a) = \lambda^{-1}\log\bigl(q_{aa}/p_{a}^{2}\bigr)$ . O triptofano é raro ( $p_{W} \approx 0.013$ ), de modo que a chance de dois triptofanos se alinharem ao acaso, $p_{W}^{2}$ , é ínfima, e um par de triptofanos conservado é um sinal de homologia muito mais forte que um par de leucinas conservado ( $p_{L} \approx 0.1$ ); a razão de log-chances é correspondentemente maior.

Exercício 5.4 ★

O que é um valor E? Uma busca devolve um acerto com $E = 3$ . O que esse número significa, e o acerto é uma homóloga?

Solução

Solução de Exercício 5.4.

O valor E é o número de alinhamentos com escore ao menos tão alto que se esperaria por acaso numa busca desta consulta contra um banco deste tamanho. $E = 3$ significa que três escores desses são esperados por acaso: o acerto não é evidência de homologia (ele ainda pode sê-lo, mas a busca não tem como dizer).

Exercício 5.5 ★★

Uma consulta de $400$ resíduos é buscada contra $2\times 10^{11}$ resíduos. Calcule o valor E de acertos com escores em bits $45$ , $55$ e $65$ . Que escore em bits dá $E = 10^{-3}$ ? Como muda a resposta se a consulta tiver $40$ resíduos?

Solução

Solução de Exercício 5.5.

$mn = 400\times 2\times 10^{11} = 8\times 10^{13} = 2^{46.2}$ . $E(45) = 2^{1.2} \approx 2.3$ ; $E(55) = 2^{-8.8} \approx 2\times 10^{-3}$ ; $E(65) = 2^{-18.8} \approx 2\times 10^{-6}$ . $E = 10^{-3}$ exige $S' = 46.2 + 10.0 = 56$ bits. Uma consulta de $40$ resíduos tem um $mn$ dez vezes menor, $2^{42.9}$ : $53$ bits bastam — mas uma consulta curta raramente chega sequer a isso.

Exercício 5.6 ★★

Calcule o conteúdo de informação de um motivo cujas quatro posições têm frequências de base (A, C, G, T) iguais a $(1,0,0,0)$ , $(0.5,0,0.5,0)$ , $(0.25,0.25,0.25,0.25)$ e $(0.7,0.1,0.1,0.1)$ . Quantas correspondências ao acaso ele tem num genoma de $4.6\,\mathrm{Mb}$ ?

Solução

Solução de Exercício 5.6.

Conteúdos de informação: $2$ , $1$ , $0$ e $2 - H$ com $H = -(0.7\log_{2} 0.7 + 3\times 0.1\log_{2} 0.1) = 0.36 + 1.00 = 1.36$ , logo $0.64$ . Total $R = 3.64$ bits. Correspondências ao acaso: $9.2\times 10^{6}$ posições em duas fitas $\times 2^{-3.64} \approx 7\times 10^{5}$ — o motivo é quase inútil sozinho.

Exercício 5.7 ★★

Explique por que as penalidades afins de lacuna são mais realistas que as lineares, e por que uma penalidade de abertura muito alta e uma muito baixa dão, as duas, alinhamentos ruins.

Solução

Solução de Exercício 5.7.

Uma inserção de vários resíduos é um único evento mutacional, de modo que seu custo não deveria crescer linearmente com seu comprimento: um custo de abertura mais um pequeno custo de extensão modela isso. Uma penalidade de abertura alta demais força discordâncias onde caberia uma lacuna e desalinha tudo depois de uma inserção verdadeira; uma penalidade baixa demais espalha lacunas por toda parte, emparelhando resíduos ao acaso e inflando a identidade.

Exercício 5.8 ★★

Uma busca BLAST de uma proteína humana contra um banco de mosca dá um melhor acerto com $E = 10^{-30}$ cobrindo os resíduos 50–180 da consulta de $600$ resíduos. A proteína da mosca é a ortóloga da humana? Que teste adicional você faria?

Solução

Solução de Exercício 5.8.

Não necessariamente: o alinhamento cobre um segmento de $130$ resíduos, que é a assinatura de um domínio compartilhado, e não de uma ortóloga alinhada em todo o seu comprimento. Teste: busque a proteína da mosca de volta contra o proteoma humano (a consulta é seu melhor acerto, em todo o comprimento?), identifique o domínio com um HMM de perfil e construa uma árvore de genes da família em várias espécies.

Exercício 5.9 ★★

Por que os métodos de perfil detectam homólogas que o alinhamento par a par não detecta? Dê um exemplo de padrão de coluna que um perfil capta e que uma sequência isolada não capta.

Solução

Solução de Exercício 5.9.

Um perfil registra, coluna por coluna, o que a família tolera: uma posição que é sempre hidrofóbica mas nunca o mesmo resíduo, um resíduo catalítico invariante, uma posição que é sempre uma lacuna em metade da família. Um alinhamento par a par pontua cada resíduo contra um único outro resíduo e não tem como saber que uma valina na posição 40 é “tão boa quanto” a isoleucina que ali está na consulta. O perfil também pondera as colunas conservadas, de modo que uma similaridade fraca concentrada onde a família é conservada se torna significativa.

Exercício 5.10 ★★★

Mostre que, sob um esquema de pontuação de escore esperado positivo, o alinhamento local de Smith–Waterman de duas sequências aleatórias longas tem um escore que cresce linearmente com o comprimento delas, e explique por que isso faz a teoria do valor E falhar. O que isso implica para alinhar DNA com correspondência $+1$ e discordância $-1$ num conteúdo GC de $60\,\%$ ?

Solução

Solução de Exercício 5.10.

Com escore esperado positivo $\mu > 0$ por coluna, o escore acumulado ao longo da diagonal de duas sequências aleatórias é um passeio aleatório com deriva positiva: após $n$ colunas ele vale cerca de $\mu n$ , de modo que o melhor alinhamento local é essencialmente a coisa toda e seu escore cresce como $\mu n$ , e não como $\log n$ . A teoria de Karlin–Altschul, que exige deriva negativa para que escores altos sejam excursões raras, não se aplica e nenhum $\lambda$ existe. Para DNA a $60\,\%$ de GC, a chance de correspondência é $2(0.3^{2}) + 2(0.2^{2}) = 0.26$ , de modo que o escore esperado é $0.26 - 0.74 = -0.48$ : ainda negativo, e a estatística vale; mas um esquema como correspondência $+1$ , discordância $-0.3$ teria esperança $+0.04$ e reportaria o genoma inteiro como um só alinhamento.

Exercício 5.11 ★★★

A tabela de Needleman–Wunsch precisa de $mn$ células de memória; para dois cromossomos de $100\,\mathrm{Mb}$ isso dá $10^{16}$ . Descreva duas ideias com que os alinhadores de genomas a evitam (sementes e encadeamento; faixas) e o que cada uma abre mão.

Solução

Solução de Exercício 5.11.

Sementes e encadeamento: encontre correspondências exatas ou quase exatas de $k$ -meros entre as duas sequências com uma tabela de dispersão, guarde as que se alinham em diagonais coerentes, encadeie-as e rode a programação dinâmica só nas lacunas entre as sementes encadeadas; abre mão de alinhamentos em regiões sem semente (trechos muito divergentes). Faixas: se as duas sequências são sabidamente quase colineares, calcule apenas as células dentro de uma faixa de largura $w$ em torno da diagonal, a um custo $wn$ em vez de $mn$ ; abre mão de qualquer alinhamento com uma inserção maior que a faixa.

Exercício 5.12 ★★★

Um HMM de busca de genes em bactérias tem estados para as três posições do códon e para o DNA não codificante. Explique como o modelo consegue distinguir sequência codificadora de não codificante sem nenhuma informação sobre códons de parada (considere o uso de códons), e por que a mesma abordagem é muito mais difícil num genoma humano.

Solução

Solução de Exercício 5.12.

A sequência codificadora tem período três: as três posições do códon têm composições de bases diferentes (a terceira é a mais enviesada), e o uso de códons é desigual em cada espécie. Um modelo com três estados codificadores em sequência, cada um emitindo bases com a composição daquela posição do códon, atribui ao DNA codificante uma probabilidade maior do que o estado não codificante atribui, ao longo de uma janela de algumas dezenas de códons, mesmo sem os códons de parada. Num genoma humano os éxons são curtos ( $150\,\mathrm{bp}$ ), separados por íntrons de quilobases, de modo que o sinal codificante é breve e interrompido; o modelo também precisa reconhecer os sítios de splicing, que são sinais fracos, e a enorme quantidade de sequência não codificante produz muitos segmentos codificantes falsos.

5.7 Problema: uma sequência do mar profundo

Problema 5.1

Problema de fim de semana — uma proteína desconhecida alinhada à mão, buscada nos bancos com sua significância calculada, seu motivo regulatório pesado em bits e seu gene conferido contra a estatística das fases de leitura aberta aleatórias, terminando no valor E do melhor acerto, nos bits de que um sítio precisa e no comprimento que uma fase de leitura precisa ter para ser acreditada

Dados: uma proteína de $300$ resíduos de um anelídeo do mar profundo. Banco de proteínas: $1.2\times 10^{11}$ resíduos. Genoma do verme: $1.6\,\mathrm{Gb}$ , $38\,\%$ GC. Pontuação dos alinhamentos à mão: correspondência $+1$ , discordância $-1$ , lacuna $-1$ . Escore em bits do melhor acerto do BLAST: $92$ ; do décimo acerto: $38$ .

Parte I — À mão.

Alinhe os peptídeos KQT e KAQT com a recorrência de Needleman–Wunsch: escreva a tabela e dê o alinhamento ótimo e o escore.
Repita com Smith–Waterman (local) para GATCAT contra ACAT: encontre o melhor alinhamento local e seu escore.
Quantas atualizações de célula custa um alinhamento global da proteína de $300$ resíduos contra uma proteína de $450$ resíduos? E contra o banco inteiro?
Se um computador faz $10^{9}$ atualizações por segundo, quanto tempo leva o alinhamento da questão 3 contra o banco inteiro? Por que se usa o BLAST em vez dele?
Um escore de identidade da BLOSUM62 é $+4$ para a alanina ( $p_{A} = 0.074$ ) e $+11$ para o triptofano ( $p_{W} = 0.013$ ). Com $\lambda = 0.347$ (unidades de meio bit), calcule a frequência-alvo $q_{AA}$ e $q_{WW}$ , e a razão $q/p^{2}$ de cada uma. Interprete.
Duas proteínas compartilham $24\,\%$ de identidade ao longo de $250$ resíduos. Diga por que a identidade sozinha não resolve a homologia aqui e o que resolveria.

Parte II — A busca.

Calcule $mn$ para a consulta contra o banco, e $\log_{2}(mn)$ .
Calcule o valor E do melhor acerto ( $S' = 92$ ) e do décimo acerto ( $S' = 38$ ).
Que escore em bits corresponde a $E = 10^{-3}$ nesta busca? E a $E = 1$ ?
O mesmo melhor acerto é encontrado quando o banco cresceu para $1.2\times 10^{12}$ resíduos. Seu valor E?
O décimo acerto alinha os resíduos 200–260 da consulta com $40\,\%$ de identidade ao longo de $60$ resíduos. Usando o valor E, diga se ele é evidência de homologia e o que uma busca por perfil poderia acrescentar.
O melhor acerto é uma cinase humana, alinhada ao longo dos resíduos 10–290. Seu melhor acerto no proteoma do verme é a consulta. O que esse teste recíproco estabelece, e o que não estabelece?

Parte III — Um motivo.

A montante do gene há um candidato a sítio de fator de transcrição de oito posições com conteúdos de informação $2, 2, 1.6, 2, 0.8, 1.2, 0.4, 0.3$ bits. Qual é o $R$ total?
Quantas correspondências ao acaso o motivo tem no genoma de $1.6\,\mathrm{Gb}$ (as duas fitas, $3.2\times 10^{9}$ posições)?
O fator regula cerca de $200$ genes. De quantos bits um motivo precisaria para especificar sozinho $200$ sítios neste genoma?
Quanto dessa falta um segundo motivo adjacente de $8$ bits poderia suprir, se os dois precisarem ocorrer juntos dentro de um espaçamento fixo?
Uma posição com frequências $(0.5, 0.5, 0, 0)$ para (A, C, G, T): calcule sua entropia e seu conteúdo de informação.
Explique, com o argumento de informação, por que os fatores de transcrição bacterianos costumam ter sítios mais longos e mais conservados que os eucariontes.

Parte IV — O próprio gene.

Em DNA aleatório de composição de bases uniforme, qual é a probabilidade de um códon ser de parada? Qual é o número esperado de códons antes que apareça um de parada (uma distribuição geométrica)?
O genoma do verme tem $38\,\%$ GC. Recalcule a probabilidade de um códon aleatório ser de parada (TAA, TAG, TGA) com as frequências de base reais, e o comprimento esperado da fase de leitura. Em que direção um conteúdo GC baixo empurra a busca de genes?
Qual é a probabilidade de uma fase de leitura aberta aleatória ter ao menos $100$ códons? E ao menos $300$ ?
No genoma de $1.6\,\mathrm{Gb}$ , seis fases em duas fitas dão cerca de $3.2\times 10^{9}$ inícios de códon. Quantas fases de leitura aberta aleatórias de ao menos $100$ códons se esperam? E de ao menos $300$ ?
Explique por que “fase de leitura aberta com mais de $100$ códons” é um localizador de genes utilizável numa bactéria, mas não neste genoma, e o que um localizador de genes eucarionte usa em vez disso.
O gene do verme tem seis éxons de $150\,\mathrm{bp}$ em média. Explique como as leituras de sequenciamento de RNA resolvem a estrutura de éxons que a sequência genômica sozinha deixa ambígua.
Resuma: o valor E do melhor acerto (questão 8), os bits necessários para especificar $200$ sítios (questão 15) e o número esperado de fases de leitura aleatórias de $300$ códons no genoma (questão 22).

Solução

Solução de Problema 5.1.

1. Linhas K, Q, T; colunas K, A, Q, T; bordas $0,-1,-2,-3,-4$ e $0,-1,-2,-3$ . Linha K: $1, 0, -1, -2$ ; linha Q: $0, 0, 1, 0$ ; linha T: $-1, -1, 0, 2$ . Ótimo $2$ : K-QT sobre KAQT. 2. Melhor escore local $3$ : CAT contra CAT (resíduos 4–6 de GATCAT com 2–4 de ACAT); ATCAT contra A-CAT também pontua $4 - 1 = 3$ . 3. $300\times 450 = 1.35\times 10^{5}$ atualizações; contra o banco, $300\times 1.2\times 10^{11} = 3.6\times 10^{13}$ . 4. $3.6\times 10^{4}$ s, dez horas por consulta; as sementes do BLAST saltam quase toda a tabela e respondem em segundos. 5. $q_{ab} = p_{a}p_{b}e^{\lambda s}$ . Alanina: $e^{1.39} = 4.0$ , $q_{AA} = 0.074^{2}\times 4.0 = 0.022$ , razão $4$ . Triptofano: $e^{3.82} = 45$ , $q_{WW} = 0.013^{2}\times 45 = 0.0077$ , razão $45$ . Um par de triptofanos alinhado é $45$ vezes mais frequente em homólogas do que por acaso, e um par de alaninas apenas quatro vezes; ainda assim os pares de alanina são mais comuns em termos absolutos, porque a alanina é comum. 6. $24\,\%$ está na zona crepuscular, em que os alinhamentos aleatórios chegam a $15\text{ a }20\,\%$ ; o valor E do alinhamento, motivos conservados nas posições certas, uma correspondência a um HMM de perfil de família conhecida ou um enovelamento compartilhado resolveriam a questão. 7. $mn = 300\times 1.2\times 10^{11} = 3.6\times 10^{13}$ ; $\log_{2}(mn) = 45.0$ . 8. $E(92) = 2^{45 - 92} = 2^{-47} \approx 7\times 10^{-15}$ ; $E(38) = 2^{7} = 128$ . 9. $E = 10^{-3}$ em $S' = 45 + 10 = 55$ bits; $E = 1$ em $45$ bits. 10. Dez vezes o $mn$ : $E \approx 7\times 10^{-14}$ , ainda esmagador. 11. Com $E = 128$ , o décimo acerto é o que o acaso produz; uma identidade de $40\,\%$ ao longo de $60$ resíduos não é evidência. Uma busca por perfil dos resíduos 200–260 contra o banco de domínios poderia mostrar se aquele segmento é um domínio conhecido, com uma estatística que falta à comparação par a par. 12. Melhores acertos recíprocos em todo o comprimento são compatíveis com uma ortologia um a um; não a provam — uma duplicação numa das linhagens depois da separação dá dois co-ortólogos, e a perda da ortóloga verdadeira pode deixar uma paráloga como melhor acerto. Uma árvore de genes com várias espécies é o teste. 13. $R = 2 + 2 + 1.6 + 2 + 0.8 + 1.2 + 0.4 + 0.3 = 10.3$ bits. 14. $3.2\times 10^{9}\times 2^{-10.3} \approx 2.5\times 10^{6}$ correspondências ao acaso. 15. $\log_{2}(3.2\times 10^{9}/200) = \log_{2}(1.6\times 10^{7}) \approx 24$ bits. 16. A coocorrência a espaçamento fixo soma os bits: $10.3 + 8 = 18.3$ , o que supre $8$ dos $13.7$ que faltam; cerca de $5.7$ bits (um fator de $50$ nas correspondências ao acaso) precisam vir de outro lugar — acessibilidade da cromatina, outros parceiros. 17. $H = -(0.5\log_{2}0.5 + 0.5\log_{2}0.5) = 1$ bit; $R = 2 - 1 = 1$ bit. 18. Um fator bacteriano precisa achar seus poucos sítios num genoma de $4.6\,\mathrm{Mb}$ sem ajuda alguma da cromatina: ele precisa de uns $19$ bits, e seus sítios são longos e conservados. Um genoma eucarionte é mil vezes maior e exigiria dez bits a mais, e ainda assim seus fatores têm sítios curtos; eles alcançam especificidade por combinação e pela restrição da cromatina acessível, o que também torna a regulação mais evoluível, já que um sítio curto é facilmente ganho ou perdido. 19. $3/64 = 0.047$ ; o número esperado de códons antes de um de parada é $64/3 \approx 21$ . 20. $p_{A} = p_{T} = 0.31$ , $p_{G} = p_{C} = 0.19$ : $P(\text{TAA}) = 0.31^{3} = 0.030$ , $P(\text{TAG}) = P(\text{TGA}) = 0.31^{2}\times 0.19 = 0.018$ ; total $0.066$ , comprimento esperado de fase $15$ códons. O DNA rico em AT é cheio de códons de parada, de modo que as fases abertas aleatórias são mais curtas e as longas se destacam mais. 21. $(61/64)^{100} = e^{-4.80} = 0.008$ ; $(61/64)^{300} = e^{-14.4} = 5.6\times 10^{-7}$ . 22. Cada fase aberta maximal termina num códon de parada, e $3.2\times 10^{9}$ inícios de códon contêm $3.2\times 10^{9}\times 3/64 = 1.5\times 10^{8}$ códons de parada: cerca de $1.5\times 10^{8}\times 0.008 = 1.2\times 10^{6}$ fases aleatórias de ao menos $100$ códons, e $1.5\times 10^{8}\times 5.6\times 10^{-7} \approx 80$ de ao menos $300$ . 23. Uma bactéria de $4.6\,\mathrm{Mb}$ tem uns $4\times 10^{5}$ códons de parada e, portanto, cerca de $3500$ fases de acaso de $100$ códons, mas quase nenhuma de $300$ ; seus genes têm em média $300$ códons e $88\,\%$ do DNA é codificante, de modo que uma fase aberta longa é quase sempre um gene. No verme, $1.5\,\%$ do DNA codifica, os éxons têm em média $50$ códons — menos que o limiar do acaso — e um milhão de fases aleatórias de $100$ códons os soterra. Os localizadores de genes eucariontes usam sinais de sítio de splicing, viés de códons num modelo oculto de Markov, homologia com proteínas conhecidas e, acima de tudo, transcritos sequenciados. 24. Uma leitura de um mensageiro processado se alinha ao genoma em dois pedaços separados por um íntron: a divisão marca os dois sítios de splicing base a base; a cobertura de leituras delineia os éxons e as leituras pareadas ligam éxons sucessivos num mesmo transcrito, resolvendo qual entre vários sítios candidatos de splicing é usado. 25. $E \approx 7\times 10^{-15}$ para o melhor acerto; cerca de $24$ bits para especificar $200$ sítios no genoma; umas $80$ fases de leitura de acaso de $300$ códons no genoma inteiro.

Termos definidos neste capítulo

Ver todos os 479 termos do glossário