Online Bioinfo
Algoritmo de Needleman-Wunsch
by Raquel Minardi in Algoritmos 18 de setembro de 2019
O algoritmo Needleman-Wunsch [1] foi proposto por Saul Needleman e Christian Wunsch na década de 1970. Trata-se de um algoritmo para alinhamento global par-a-par de sequências. Isso significa que ele alinha pares de sequências em todo o seu comprimento. É importante ainda destacar que ele é ótimo, ou seja, retorna o alinhamento de maior pontuação de acordo com o esquema de pontuação utilizado. É interessante notar que esse algoritmo foi proposto antes da publicação das matrizes de substituição PAM e BLOSUM e utiliza um esquema de pontuação genérico baseado na distância de edição de sequências.
O algoritmo de Needleman-Wunsch foi descoberto e redescoberto inúmeras vezes em áreas que vão desde reconhecimento de fala até, obviamente, biologia molecular [2]. Apesar de haverem detalhes que diferem nas diversas aplicações, todas são essencialmente programação dinâmica.
Problema do Turista em Manhattan
Antes de entrar no alinhamento de sequências, mostraremos um outro problema que também é resolvido da mesma forma. Esse problema se chama problema do Turista em Manhattan e servirá para desenvolver sua intuição para que compreenda mais facilmente o problema do alinhamento de sequências.
Apresentamos o mapa do centro de Manhattan onde se pode ver inúmeros pontos turísticos como o Carnegie Hall, o Rockfeller Center, a Times Square, entre outros. Uma característica interessante dessa região que se pode notar através do mapa é a organização bem retangular dos quarteirões. Imagine que você está em Manhattan e que está em uma excursão que lhe deixará na esquina entre a 59th Street e a 8th Avenue e lhe buscará no Chrysler Building na esquina entre a 42nd Street e a Lexington Avenue. Há inúmeras atracões turísticas na região entre esses pontos e você deseja organizar seu trajeto de forma a visitar o máximo de pontos turísticos quanto possível sem se preocupar com o tempo.
Como todos os quarteirões são retangulares, essa região pode ser facilmente representada como um grid como o da Figura 1 a seguir no qual temos nós (círculos) representando as esquinas e arestas (setas) representando os trechos das ruas entre as esquinas. Os rótulos das arestas indicam quantos pontos turísticas há naquele trecho da rua.
Chamaremos o primeiro nó do canto superior esquerdo de fonte e o último do canto inferior direito de sumidouro. O objetivo do turista é sair da fonte e chegar no sumidouro com o maior número possível de pontos turísticos visitados. Note que ele só pode caminhar de cima para baixo e da esquerda para a direita não sendo permitidas voltas.
Como você resolveria esse problema? Usar um método de força bruta seria uma possibilidade e consistiria em buscar todos os possíveis caminhos entre fonte e sumidouro. Essa solução seria proibitiva em termos de custo computacional, que seria fatorial. Poderíamos então usar uma abordagem gulosa que consiste em, a cada esquina, olhar as ruas adiante escolhendo a opção aparentemente mais vantajosa localmente dentre as duas possíveis ruas a tomar. Essa solução pode parecer interessante em determinado momento mas pode levar a uma solução péssima levando o turista para um local com poucas ou nenhuma atração e levando a uma solução ruim.
Felizmente, esse problema tem subestrutura ótima e pode ser resolvido por programação dinâmica. Ao invés de solucionar o problema de ir da fonte (0,0) ao sumidouro (n,m), resolveremos o problema mais geral e encontrar o caminho máximo entre a fonte e um vértice arbitrário (i,j), com 0 ≤ i ≤ n, 0 ≤ j ≤ m. Denotaremos o menor caminho por si,j, se do sn,m a solução do Problema do Turista em Manhattan.
Aparentemente, estamos escolhendo um problema geral ainda mais difícil mas isso não é verdade pois, para resolver o Problema do Turista em Manhattan, precisaremos resolver os subproblemas por programação dinâmica, então, é a mesma dificuldade.
Comecemos por s0,j (para 0 ≤ j ≤ m) que é um problema simples bastando que o turista se mova sempre à direita e acumulando as pontuações. A seguir, podemos fazer o mesmo para si,0 (para 0 ≤ i ≤ n) em que o turista se moverá sempre para baixo e acumulando as pontuações. Note que colocamos a soma da pontuação do caminho da fonte até o ponto si,j como um rótulo dentro do nó.
A seguir, podemos preencher tanto a linha 2 quando a coluna 2 e assim sucessivamente. Note que as novas soluções que vão sendo geradas para cada problema (cada nó é um tamanho de problema) são compostas pelas soluções dos subproblemas menores que já foram calculados. Por exemplo, a célula s1,1 usará as soluções s0,1 ou s1,0. Nesse caso, a melhor solução será vir de s1,0 somando o valor da aresta que é 3. Preenchemos as linhas e colunas sucessivamente.
Finalmente, termos a matriz completa com os caminhos máximos saindo da fonte s0,0 para todos os possíveis si,j com 0 ≤ i ≤ n, 0 ≤ j ≤ m teremos:
Interessantemente, essa matriz pode ser preenchida em três possíveis sequências conforme ilustrado nas Figuras 4, 5 e 6. Não importa a ordem em que seja preenchida, seu conteúdo será o mesmo.
Voltando ao alinhamento de sequências...
O alinhamento de duas sequências v (com n caracteres) e w (com m caracteres) pode ser visualizado assim:
Note que esse alinhamento poderá ter, no máximo (n+m) posições. As colunas que contêm o mesmo caracter são chamadas de matches (combinações). As colunas que contêm espaços (gaps) são chamadas indels (inserções e deleções) sendo que quando a linha superior contem um caracter é uma deleção e quando ela contêm um espaço, é uma inserção.
O alinhamento exemplificado acima tem 5 matches (verde), 0 mismatches (pareamento de caracteres não idênticos) e 4 indels (amarelo). Esse alinhamento poderia ser representado como o seguinte grid ou matriz na qual a primeira sequência está representada nas colunas e a segunda nas linhas (Figura 8):
Essa matriz é central em todo algoritmo de alinhamento e o caminho destacado na Figura 9 representa o melhor alinhamento entre s0,0 e sn,m. É interessante notar que as arestas no caminho máximo representam impressões:
Horizontais: significam deleções na primeira sequência v, ou seja, um caracter na primeira sequência é impresso e um gap é impresso na segunda sequência w.
Verticais: significam inserções na primeira sequência v, ou seja, um gap é impresso na primeira sequência e um caracter é impresso na segunda sequência w.
Diagonais: significam um match ou mismatch, ou seja, um caracter é impresso em cada uma das sequências v e w podendo ser iguais ou diferentes.
Falta então um esquema de pontuação que seja capaz de julgar o mérito dos diversos possíveis alinhamentos. Teremos, então, o custo de admitir matches e mismatches bem como dos indels.
Suponha, por simplicidade, que usaremos um esquema demasiadamente simplificado no qual pontuamos "+1" para um match e "0", caso contrário. Nessa primeira versão do algoritmo, não contemplaremos um mismatch visto que o esquema de pontuação permite caminhamento em diagonal apenas no caso de igualdade entre os caracteres. Destacamos que não há perda de generalidade pois bastaria criar uma pontuação diferenciada para mismatches e o mesmo algoritmo poderia ser aplicado.
Esse problema é clássico em ciência da computação e é chamado de Problema da Máxima Subsequência Comum (do inglês Longest Common Subsequence - LCS) e iremos abordá-lo como um aquecimento.
Problema da Máxima Subsequência Comum
Encontre a máxima subsequência comum entre duas sequências
Entradas: Duas sequências, v e w
Saída: A máxima subsequência comum entre v e w
Para resolver o LCS, criaremos uma matriz de programação dinâmica colocando a primeira sequência como as colunas e a segunda como as linhas conforme a seguir e seguindo o seguinte critério si, j = max (si-1,j, si,j-1 e si-1,j-1+1 (desde que v[i] = w[j]):
Note que preenchemos a primeira linha (s0,j) e coluna (si,0) com 0s para inicializar a pontuação do alinhamento. Vamos preenchendo essas matrizes com as pontuações obtidas nos possíveis alinhamentos. Comecemos com a célula s1,1 que corresponde ao problema de se alinhar as sequências "A" com "A". Nesse caso, teríamos um match e pontuaríamos "+1" e deixando uma marca diagonal.
Podemos prosseguir, calculando os problemas de toda a linha s1,j. Note que a célula s1,2, representaria um mismatch de "A" com "T" e pontuaria "0". Essa pontuação tem de ser somada com a de s1,1 pois a posição corrente indica o alinhamento de "A" com "AT" e valeria "1" com a marcação de horizontal pois veio de s1,1. Por fim, podemos continuar o preenchimento da matriz até o fim:
Dessa matriz, podemos perceber que o alinhamento de pontuação máxima terá valor 5. Além disso, seguindo os ponteiros partindo da célula sn,m até chegar a célula s0,0, podemos gerar o alinhamento:
Note que esse alinhamento é construído de trás para frente, ou seja, da direita para a esquerda. Partindo da célula s7,7, na qual temos um "5←", imprimiremos no alinhamento o caracter "C" da sequência v e um gap na sequência w. Recapitulando, essa decisão se deve à direção apontada pelo ponteiro:
Horizontais (←): significam deleções na primeira sequência v, ou seja, um caracter na primeira sequência é impresso e um gap é impresso na segunda sequência w.
Verticais (↑): significam inserções na primeira sequência v, ou seja, um gap é impresso na primeira sequência e um caracter é impresso na segunda sequência w.
Diagonais (↖): significam um match, ou seja, um caracter é impresso em cada uma das sequências v e w.
Por fim, lembramos que o algoritmo de Needleman-Wunsch realiza um alinhamento global par-a-par de sequências que significa que ele alinha pares de sequências em todo o seu comprimento. Seu resultado é ótimo, ou seja, retorna o alinhamento de maior pontuação de acordo com o esquema de pontuação utilizado. Contudo, é importante destacar que pode existir mais de um alinhamento com a mesma pontuação e, nesse caso, um deles será retornado. A escolha de qual alinhamento é gerado depende da prioridade dada às condições do condicional para preenchimento da matriz de programação dinâmica da implementação do algoritmo.
[1] Needleman, S.B. e Wunsch, C.D.. "A general method applicable to the search for similarities in the amino acid sequence of two proteins." Journal of molecular biology 48.3 (1970): 443-453.