SVD – Decomposição por Valores Singulares com Python

Este material faz parte do livro “Python para Bioinformática” (ainda não publicado). Para mais informações, consulte o autor.

Decomposição por valores singulares, do inglês Singular Value Decomposition (SVD) é uma técnica da álgebra linear usada para recuperação de informações, redução dimensional e remoção de ruídos. Nessa técnica, uma matriz A de tamanho MxN é decomposta em três outras matrizes, aqui denominadas como U, S e V.

A=U*S*V’

Sendo as U e V matrizes ortonormais, de dimensões MxM e NxN, respectivamente; e S a matriz diagonal com os valores singulares.

Neste livro, não abordaremos a fundo os pilares matemáticos que descrevem a teoria de como essa técnica funciona. Por isso, tenha em vista apenas que, em machine learning e bioinformática, a técnica de SVD tem diversas aplicações. Na seção a seguir, faremos uma comparação entre sequências usando matrizes esparsas e SVD.

Comparando sequências usando SVD

Sei que já falamos muito neste livro sobre comparações de sequências. Logo, você deve estar se perguntando por que precisa aprender um novo método? Além disso, o método que usaremos nesta seção não é melhor que os métodos tradicionais baseados em alinhamentos. O que faremos é tentar mostrar que eles podem apresentar um resultado levemente próximo.

Então por que usar o SVD? De fato, você talvez não precise aprender um novo método de comparação entre sequências. Todavia, tenha em vista que a comparação que faremos aqui será apenas uma desculpa para que possamos compreender e aplicar os fundamentos do SVD.

Em certos contextos, o uso do SVD pode acelerar muito sua aplicação, uma vez que seu uso permite a redução dimensional de grandes volumes de dados. Entretanto, aplicar isso na prática é bastante complexo. O que faremos aqui é tentar exemplificar um contexto de aplicação dessa técnica. Assim, conforme você for compreendendo os conceitos, poderá aplicá-los em outros contextos.

Passo 1 – converta a sequência em uma matriz binária esparsa

Para que possamos usar a técnica de SVD, precisamos ter como entrada uma matriz binária esparsa. De uma maneira bem leiga, podemos dizer que uma matriz esparsa é uma tabela composta por números 0 ou 1. Então, como podemos converter uma sequência em uma matriz de 0s ou 1s?

Para isso, vamos usar uma variação da técnica de janela deslizante (K-mer), que apresentamos anteriormente. Na técnica de janela deslizante, dividimos a sequência em pequenos fragmentos de tamanho K, que se sobrepõe. Por exemplo, a sequência de nucleotídeos ATCAT poderia ser quebrada em três sequências de tamanho K=3 (ATC, TCA e CAT). Veja:

ATCAT
ATC  
 TCA
  CAT

Agora se quisermos quebrar nossa sequência em subsequências de tamanho K=2, teríamos quatro subsequências (AT, TC, CA, AT):

ATCAT
AT  
 TC
  CA
   AT

Chamamos cada subsequência de k-mer (ou 2-mer, neste caso).

Agora, se quisermos converter essa sequência em uma matriz esparsa, precisaremos criar um campo para cada possível K-mer (vamos usar o tamanho K=2 como exemplo) e, em seguida, verificar se ele existe na sequência.

Para K=2, temos 16 possíveis combinações (ou seja, AA, AC, […], TT). Assim, devemos criar uma tabela com as combinações, além de uma coluna para dizer se o k-mer existe na sequência em questão. Se existir, vamos preencher esse campo com 1. Se não existir, vamos preencher esse campo com 0. Observe como ficará essa tabela para a sequência ATCAT com K=2:

#	Combinações (K=2)	Existe em ATCAT
1	AA	0
2	AC	0
3	AG	0
4	AT	1
5	CA	1
6	CC	0
7	CG	0
8	CT	0
9	GA	0
10	GC	0
11	GG	0
12	GT	0
13	TA	0
14	TC	1
15	TG	0
16	TT	0

Viu como podemos representar uma sequência como uma matriz binária esparsa? Neste caso, apenas três campos foram preenchidos com 1 (note que o k-mer “AT” aparece duas vezes, mas ignoramos isso). O resto dos campos foram preenchidos com 0. Simples, não é?

Agora, se quisermos fazer isso com uma sequência de proteínas será um pouco mais complicado, pois temos 20 aminoácidos. Por exemplo, se quisermos analisar um K=2, teríamos 400 combinações. Para K=3, teríamos 8000 combinações.

Vamos construir um programa que gera a matriz esparsa binária.

Primeiro vamos receber os dados de entrada:

# dados de entrada
entrada = "ATCAT"
k = 2
n = len(entrada)

Note que, antes de calcular os k-mers, precisamos definir o tamanho de k e precisamos ainda obter o tamanho da sequência.

Agora, vamos criar uma função que determina todas as combinações possíveis. Atenção: para isso, precisaremos criar uma função recursiva. Recursividade é um assunto que pode nos dar calafrios. Isso porque uma função recursiva precisa chamar a si mesmo várias e várias vezes. Observe como ficará o código e depois vamos discuti-lo:

 5. # PARTE 1 - armazena todas as k combinações
 6. def combinador(k, letras):
 7.   if k == 0:
 8.     return ['']
 9.     
10.     combinacoes = []
11.     combinacoes_anteriores = combinador(k - 1, letras)
12.       
13.     for letra in letras:
14.       for combinacao_anterior in combinacoes_anteriores:
15.         combinacoes.append(letra + combinacao_anterior)
16.    
17.     return combinacoes
18.    
19.   combinacoes = combinador(k, "ATCG")
20.   # print(combinacoes)

Nossa função “combinador” recebe como entrada: o valor de k (neste exemplo, 2) e todas as letras que serão usadas nas combinações (primeiramente, vamos fazer para os nucleotídeos “ATCG” e depois refaremos o exemplo para aminoácidos). Observe que no início da função, verificamos se k é igual a zero e retornamos um vetor com uma string vazia. Esse trecho é a condição de parada do nosso algoritmo recursivo (só por curiosidade, em computação chamamos isso de “caso base”). Note que na linha 10 criamos uma lista para armazenar nossas combinações (isso também servirá para resetar o valor armazenado nessa lista). Em seguida, na linha 11, criamos uma lista para armazenar as combinações anteriores.

Mas espere!

A lista de combinações anteriores chama a própria função “combinador”. É nesse ponto que nosso código começa a ficar confuso. Mas não se assuste. Funções recursivas precisam chamar elas mesmas passando como entrada um fragmento menor do problema. Nesse caso, estamos fazendo uma chamada enviando um valor menor de k.

Agora, você entende por que precisamos de uma condição de parada que retorna um vetor com uma string vazia? A cada vez que uma função recursiva chama a si mesmo com outros valores de entrada, essa nova execução tem prioridade e a execução atual entra em uma espécie de “fila de espera” para ser concluída.

Para ilustrar isso, vamos dizer que chamamos o combinador para k=3.

combinador(3, "ATCG")

logo, o seguinte código será executado (k=3):

if k == 0:
  return ['']
combinacoes = []
combinacoes_anteriores = combinador(2, letras) # entrada: k-1

Neste ponto, o próprio código executa o combinador de novo, passando como entrada o valor 2. Assim, o seguinte código será executado nessa nova chamada (k=2):

if k == 0:
  return ['']
combinacoes = []
combinacoes_anteriores = combinador(1, letras) # entrada: k-1

Mais uma vez o combinador será executado, agora passando como entrada o valor 1 (k=1):

if k == 0:
  return ['']
combinacoes = []
combinacoes_anteriores = combinador(0, letras) # entrada: k-1

Ao chegar nessa linha, o combinador é executado novamente passando como entrada o valor 0:

if k == 0:
  return ['']

Espere! Agora k vale zero, então ele retorna o seguinte valor: [“”].

Agora sim ele continuará a execução. Mas preste atenção que há várias funções esperando para serem executadas. Elas serão chamadas de acordo com sua posição na fila:

·   combinacoes (k=1)

·   combinacoes (k=2)

·   combinacoes (k=3)

Nesse exemplo, a próxima execução da função de combinações ocorre quando k=1. Preste atenção nos valores presentes na variável “combinacoes_anteriores”. Esses valores são a chave para entender como o algoritmo recursivo funciona (para facilitar, vamos repetir o código adicionando o valor dessa variável em um comentário logo a frente).

for letra in letras: #["A","T","C","G"]
  for combinacao_anterior in combinacoes_anteriores: #[""]
    combinacoes.append(letra + combinacao_anterior)
return combinacoes

Nesse momento, letras tem quatro itens: “ATCG”. Enquanto isso, combinacoes_anteriores possui apenas uma lista vazia [‘’].

Neste caso, a variável “combinacoes” vai receber os seguintes valores:

[“A”, “T”, “C”, “G”]

Note que usamos return para sair da execução da função. Agora precisamos continuar executando as funções que estavam na fila de espera:

·  combinacoes (k=1)

·  combinacoes (k=2)

·  combinacoes (k=3)

Ou seja, será executado o código para k=2:

for letra in letras: # ["A","T","C","G"]
  for combinacao_anterior in combinacoes_anteriores: # ["A","T","C","G"]
    combinacoes.append(letra + combinacao_anterior)
return combinacoes

Agora, a lista de combinações anteriores tem quatro valores [“A”,”T”,”C”,”G”]. Logo, a lista de combinações (que foi zerada anteriormente) receberá os seguintes valores:

['AA', 'AT', 'AC', 'AG', 'TA', 'TT', 'TC', 'TG', 'CA', 'CT', 'CC', 'CG', 'GA', 'GT', 'GC', 'GG']

Por fim, vamos para a última função que estava aguardando na fila:

·  combinacoes (k=1)

·  combinacoes (k=2)

·  combinacoes (k=3)

for letra in letras: # ["A","T","C","G"]
  for combinacao_anterior in combinacoes_anteriores: # ['AA', 'AT', 'AC', 'AG', 'TA', 'TT', 'TC', 'TG', 'CA', 'CT', 'CC', 'CG', 'GA', 'GT', 'GC', 'GG']
    combinacoes.append(letra + combinacao_anterior)
return combinacoes

Perceba que a lista de combinações anteriores possui 16 valores. Quando ocorrer a combinação com os quatro valores de letras, teremos 64 combinações de tamanho 3:

['AAA', 'AAT', 'AAC', 'AAG', 'ATA', 'ATT', 'ATC', 'ATG', 'ACA', 'ACT', 'ACC', 'ACG', 'AGA', 'AGT', 'AGC', 'AGG', 'TAA', 'TAT', 'TAC', 'TAG', 'TTA', 'TTT', 'TTC', 'TTG', 'TCA', 'TCT', 'TCC', 'TCG', 'TGA', 'TGT', 'TGC', 'TGG', 'CAA', 'CAT', 'CAC', 'CAG', 'CTA', 'CTT', 'CTC', 'CTG', 'CCA', 'CCT', 'CCC', 'CCG', 'CGA', 'CGT', 'CGC', 'CGG', 'GAA', 'GAT', 'GAC', 'GAG', 'GTA', 'GTT', 'GTC', 'GTG', 'GCA', 'GCT', 'GCC', 'GCG', 'GGA', 'GGT', 'GGC', 'GGG']

Agora que todas as funções recursivas foram executadas e a variável “combinacoes” armzena os k-mers existentes, precisamos criar a matriz esparsa e preenchê-la com zeros. Podemos fazer isso usando um dicionário. Veja:

21.   # PARTE 2 - cria uma matriz vazia
22.   matriz = {}
23.   for i in combinacoes:
24.     matriz[i] = 0
25.   # print(matriz)

Veja um exemplo da saída para k=2:

{'AA': 0, 'AT': 0, 'AC': 0, 'AG': 0, 'TA': 0, 'TT': 0, 'TC': 0, 'TG': 0, 'CA': 0, 'CT': 0, 'CC': 0, 'CG': 0, 'GA': 0, 'GT': 0, 'GC': 0, 'GG': 0}

Agora sim podemos ler a sequência de entrada e verificar cada k-mer dela. Lembre-se que a quantidade de k-mers será sempre o tamanho da sequência menos o tamanho de k mais 1. Veja como fica o código:

26.   # PARTE 3 - preenche a matriz esparsa
27.   for i in range(0, n-k+1):
28.     k_mer = entrada[i:i+k]
29.     matriz[k_mer] = 1
30.    
31.   print(matriz)

Para o exemplo apresentado, nosso código imprimiria:

{'AA': 0, 'AT': 1, 'AC': 0, 'AG': 0, 'TA': 0, 'TT': 0, 'TC': 1, 'TG': 0, 'CA': 1, 'CT': 0, 'CC': 0, 'CG': 0, 'GA': 0, 'GT': 0, 'GC': 0, 'GG': 0}

Agora vamos unir todos os códigos que construímos e testar com sequências de proteínas. Note que a única grande mudança é que nossas letras não serão mais apenas “ATCG”, mas sim as 20 letras dos aminoácidos: “ACDEFGHIKLMNPQRSTVWY”. Ou seja, basta mudar a chamada inicial do combinador:

combinacoes = combinador(k, "ACDEFGHIKLMNPQRSTVWY")

Vamos usar como entrada, as três sequências de hemoglobinas:

Humano	Boi	Chimpanzé
MVHLTPEEKSAVTALWGKVNVD EVGGEALGRLLVVYPWTQRFFES FGDLSTPDAVMGNPKVKAHGKK VLGAFSDGLAHLDNLKGTFATLS ELHCDKLHVDPENFRLLGNVLVC VLAHHFGKEFTPPVQAAYQKVV AGVANALAHKYH	MVLSAADKGNVKAAWGKVGGHA AEYGAEALERMFLSFPTTKTYFPH FDLSHGSAQVKGHGAKVAAALTK AVEHLDDLPGALSELSDLHAHKLR VDPVNFKLLSHSLLVTLASHLPSD FTPAVHASLDKFLANVSTVLTSKY R	MVHFTAEEKAAVTSLWSKM NVEEAGGEALGRLLVVYPWT QRFFDSFGNLSSPSAILGNPK VKAHGKKVLTSFGDAIKNMD NLKPAFAKLSELHCDKLHVDP ENFKLLGNVMVIILATHFGKEF TPEVQAAWQKLVSAVAIALAH KYH

Para facilitar a análise desse código, vamos dividir as principais partes do código em funções. Assim poderemos reutilizá-las várias vezes sem a necessidade de replicar várias vezes os mesmos códigos (preste atenção que, ao trabalhar com matrizes em funções, é necessário fazer uma cópia delas usando o método copy).

Observe como ficará o código completo:

 1. humano = "MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH"
 2.  
 3. boi = "MVLSAADKGNVKAAWGKVGGHAAEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHGAKVAAALTKAVEHLDDLPGALSELSDLHAHKLRVDPVNFKLLSHSLLVTLASHLPSDFTPAVHASLDKFLANVSTVLTSKYR"
 4.  
 5. chimpanze = "MVHFTAEEKAAVTSLWSKMNVEEAGGEALGRLLVVYPWTQRFFDSFGNLSSPSAILGNPKVKAHGKKVLTSFGDAIKNMDNLKPAFAKLSELHCDKLHVDPENFKLLGNVMVIILATHFGKEFTPEVQAAWQKLVSAVAIALAHKYH"
 6.  
 7. # PARTE 1- armazena todas as k combinações
 8. def combinador(k, letras):
 9.   if k == 0:
10.       return ['']
11.       
12.     combinacoes = []
13.     combinacoes_anteriores = combinador(k - 1, letras)
14.      
15.     for letra in letras:
16.       for combinacao_anterior in combinacoes_anteriores:
17.         combinacoes.append(letra + combinacao_anterior)
18.    
19.     return combinacoes
20.    
21.   # PARTE 2 - cria uma matriz vazia
22.   def cria_matriz_zeros(combinacoes):
23.     matriz = {}
24.     for i in combinacoes:
25.       matriz[i] = 0
26.    
27.     return matriz
28.    
29.   # PARTE 3 - preenche a matriz esparsa
30.   def preenche_matriz_binaria(entrada, matriz):
31.     n = len(entrada)
32.    
33.     matriz_backup = matriz.copy()
34.    
35.     for i in range(0, n-k+1):
36.       k_mer = entrada[i:i+k]
37.       matriz_backup[k_mer] = 1
38.    
39.     return matriz_backup
40.    
41.   # dados de entrada
42.   k = 2
43.   combinacoes = combinador(k, "ACDEFGHIKLMNPQRSTVWY")
44.   matriz = cria_matriz_zeros(combinacoes)
45.    
46.   # cria as matrizes
47.   matriz_humano = preenche_matriz_binaria(humano, matriz)
48.   matriz_boi = preenche_matriz_binaria(boi, matriz)
49.   matriz_chimpanze = preenche_matriz_binaria(chimpanze, matriz)
50.    
51.   # imprimindo as matrizes
52.   print(matriz_humano) # {'AA': 1, 'AC': 0, 'AD': 0, 'AE': 0, [...]
53.   print(matriz_boi)   # {'AA': 1, 'AC': 0, 'AD': 1, 'AE': 1, [...]
54.   print(matriz_chimpanze) # {'AA': 1, 'AC': 0, 'AD': 0, 'AE': 1, [...]

Note que, até o momento, não temos uma “matriz esparsa binária” de verdade. Temos apenas três dicionários com dados. Assim, vamos unir esses valores em uma única matriz (para dar um toque dramático, a chamarei de “matrix”).

55. matrix = []
56. 
57. matrix.append(list(matriz_humano.values()))
58. matrix.append(list(matriz_boi.values()))
59. matrix.append(list(matriz_chimpanze.values()))
60. 
61. print("Linhas:", len(matrix), "\nColunas:", len(matrix[0]))
62. # Linhas: 3
63. # Colunas: 400

Nossa “matrix” possui 3 linhas por 400 colunas. Não vamos imprimi-la aqui por completo, pois são muitos 1s e 0s, mas se você a imprimi-la, verá algo parecido com isto:

[
 [1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, ... 0],
 [1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, ... 0],
 [1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, ... 0]
]

Os k-mers não são tão relevantes para o que faremos a seguir, por isso os removemos chamando o método “values()” e aplicando a função “list()”.

Agora que já temos uma matriz binária para cada sequência podemos começar a usar o SVD para compará-las.

Passo 2 – Decompondo uma matriz usando SVD

Agora, vamos usar o SVD para realizar uma redução dimensional. Veja que, para esse estudo de caso, partimos de três sequências de aproximadamente 150 letras cada. Em seguida, convertemos isso em três listas com 400 números. Agora, usaremos o SVD para converter esses 400 valores que representam cada proteína em apenas dois (convenhamos que é bem mais fácil comparar dois números do que 150 letras ou 400 valores entre 0s ou 1s). O que queremos aqui é mostrar que apenas esses dois números serão suficientes para representar cada uma das três proteínas em determinado contexto.

Para usar o SVD, você precisará de três bibliotecas: numpy, matplotlib e scipy. Nos capítulos anteriores, você conheceu um pouco sobre o NumPy e sobre o Matplotlib. Scipy (https://scipy.org) é uma biblioteca Python usada para aplicações científicas que requerem análises eficientes de matrizes multidimensionais.

Para instalar a biblioteca scipy use o comando:

pip3 install scipy

Vamos começar importando as três bibliotecas:

1.  import numpy as np
2.  from scipy.linalg import svd
3.  import matplotlib.pyplot as plt

Note que usaremos apenas o método svd do scipy. Caso deseje realizar uma decomposição parcial em valor singular de uma matriz esparsa, importe a função svds do scipy.sparse.linalg.

Agora, vamos decompor nossa matriz nas três matrizes: U, S e V.

Você pode fazer isso da seguinte forma:

4.  [U, S, V] = svd(matrix)
5.  print(U.shape, S.shape, V.shape)

Observe que será impresso na tela os tamanhos das matrizes:

(3, 3) (3,) (400, 400)

Nossa matriz U tem 3 linhas por 3 colunas. A matriz S possui apenas três valores em uma única linha (note que S era para ser uma matriz diagonal, mas o scipy já extraiu os valores relevantes e os salvou em uma única linha). Por fim, V tem 400 linhas por 400 colunas.

Note que o valor 3 se refere às três sequências, enquanto 400 se refere a quantidade de K-mers existentes (como usamos k=2, temos 20² combinações).

Agora, vamos realizar uma redução dimensional. Para isso, vamos extrair duas dimensões de valores singulares e salvá-las em outras variáveis. Observe como podemos fazer isso:

6.  S_2d = S[0:2]
7.  U_2d = U[0:2, :]
8.  matrix_reduzida = S_2d*U_2d.transpose()

Observe que fizemos uso do método transpose(). Esse método transpõe uma determinada matriz, isto é, inverte colunas e linhas.

Agora que já extraímos as duas dimensões, vamos plotar na tela um ponto para representar cada sequência. Podemos fazer isso usando matplotlib.

 9.   # Criando a figura
10.   fig = plt.figure()
11.   fig.suptitle('Hemoglobinas')
12.    
13.   # extraindo os dados para os eixos x e y
14.   dados = matrix_reduzida.transpose()
15.   x = dados[0]
16.   y = dados[1]
17.    
18.   # criando um gráfico de pontos
19.   plt.scatter(x, y, s=500)
20.    
21.   # adicionando legenda
22.   legenda = ['Humano', 'Boi', "Chimpanze"]
23.   for i, especie in enumerate(legenda):
24.     plt.annotate(especie, (x[i], y[i]))
25.    
26.   # plotando na tela
27.   plt.show()

Esse código irá gerar a seguinte figura:

Aparentemente, as esferas que representam sequências de hemoglobinas de humanos e chimpanzés estão mais próximas (o que já esperávamos).

Se você quiser, podemos encontrar valores numéricos que indicam isso. Para isso, vamos construir uma função que calcula a distância euclidiana em um ambiente bidimensional. Observe como ficará o código:

28. def distancia(x1,y1,x2,y2):
29.  return ((x1-x2)**2 + (y1-y2)**2)**(1/2)
30. 
31. dist_h_b = distancia(x[0], y[0], x[1], y[1])
32. dist_h_c = distancia(x[0], y[0], x[2], y[2])
33. dist_b_c = distancia(x[1], y[1], x[2], y[2])
34. 
35. print("Humano x boi:", int(dist_h_b))   # 15
36. print("Humano x chimpanze:", int(dist_h_c)) # 3
37. print("Boi x chimpanze:", int(dist_b_c)) # 17

Observe que será impresso na tela os valores: 15, 3 e 17 (não há bem uma unidade de medida para isso). Ou seja, a distância dos pontos das sequências de hemoglobinas:

• humanas para as do boi foi de 15 (a identidade das sequências usando BLAST é de 43%);
• do humano para o chimpanzé foi de 3 (identidade das sequências é de 76%) e
• do boi para o chimpanzé foi de 17 (identidade é de 38%).

Ou seja, quanto menor a distância, maior a identidade. Isso indica que o uso do SVD como técnica de redução dimensional foi capaz de reduzir a quantidade de dados (note que agora precisamos de apenas dois números para representar cada sequência: as posições x e y), mantendo as similaridades e diferenças entre elas à vista.

Comparando com outras entradas

À medida que a quantidade de elementos aumenta, o custo para decompor a matriz usando SVD também aumenta. Por isso, decompor matrizes com muitas sequências pode levar um grande tempo. Entretanto, podemos calcular um novo elemento sem a necessidade de recalcular toda a matriz com SVD. Para isso, basta que nossa matriz original já tenha uma boa quantidade de indivíduos representativos (nesse estudo de caso infelizmente não temos, mas para fins didáticos, vamos fingir que nossas três sequências de hemoglobina representam uma boa amostra das hemoglobinas que existem pelo universo).

Vamos adicionar uma nova hemoglobina a nossa lista: a da espécie Sus scrofa, um tipo de porco (ou javali). Ela pode ser obtida no UniProt (ID: P02009).

> P02009
SLTKAERTIIGSMWTKISSQADTIGTETLERLFASYPQAKTYFPHFDLNPGSAQLRAHGSKVLAAVGEAVKSIDNVSAALAKLSELHAYVLRVDPVNFKFLSHCLLVTLASHFPADLTAEAHAAWDKFLTIVSGVLTEKYR

Vamos calcular a matriz binária de k-mer dessa sequência e então compará-la com as outras:

38. porco = "SLTKAERTIIGSMWTKISSQADTIGTETLERLFASYPQAKTYFPHFDLNPGSAQLRAHGSKVLAAVGEAVKSIDNVSAALAKLSELHAYVLRVDPVNFKFLSHCLLVTLASHFPADLTAEAHAAWDKFLTIVSGVLTEKYR"
39. matriz_porco = preenche_matriz_binaria(porco, matriz)
40. matriz_porco = list(matriz_porco.values())
41. 
42. # Extraindo 2D da matriz V
43. V_2d = V[0:2, :]
44. 
45. # Multiplicando a matriz do porco com 2D da matriz V
46. ponto_porco = np.dot(V_2d, matriz_porco)
47. 
48. # Plotando novo item
49. plt.scatter(ponto_porco[0], ponto_porco[1], s=500)
50. plt.annotate("Porco", (ponto_porco[0], ponto_porco[1]))
51. plt.show()

Observe o resultado:

Podemos ver que a sequência da hemoglobina do porco está mais próxima da do boi que a do humano e do chimpanzé estão. Podemos conferir isso checando as distâncias:

print("Porco x Boi:", int(distancia(ponto_porco[0], ponto_porco[1], x[1], y[1])))

• Porco x Boi: 11
• Humano x Boi: 15 (calculada anteriormente)
• Chimpanzé x Boi: 17 (calculada anteriormente)

O que condiz com os resultados de análises de comparações feitas com BLAST:

• Porco x Boi (56% de identidade)
• Boi x Humano (42% de identidade)

Análise de POSTO

Nos exemplos anteriores, usamos apenas duas dimensões. Entretanto, podemos usar mais dimensões caso desejarmos. Para fazer isso, podemos realizar uma análise de posto (factor). Essa análise é feita na matriz diagonal S, que é responsável por armazenar os valores singulares.

Vamos começar com uma análise do gráfico de posto:

s = S
s = s*(1/np.sum(s))
fig = plt.figure()
fig.suptitle('Gráfico de POSTO')
plt.plot(s)
plt.show()

Esse gráfico plota a curva da diagonal S. Em geral, dizemos que posto ideal é representado pelo ponto do eixo X em que se tem a maior inclinação da curva – chamamos isso de análise do cotovelo (elbow). Observe a figura plotada no exemplo:

Podemos ver que nossa curva se inclina no ponto 1.00 do eixo X do gráfico. Após isso, ela tende a entrar em um platô (dificilmente observaremos isso nessa figura, pois só temos três dimensões). Como nosso eixo X começa de 0, o ponto 1.00 indica a 2ª dimensão.

Como a análise do “cotovelo” do gráfico de posto pode ser um pouco desafiador para maioria de nós, podemos usar uma outra estratégia que usa valores numéricos. Para isso, precisamos calcular o somatório acumulativo total de valores na matriz S e determinar qual posto alcança 70% desse valor.

Observe como isso pode ser implementado:

somatorio = 0
total = sum(S)
posto_ideal_nao_descoberto = True

for i in range(len(S)):
  percentual = 100*S[i]/total
  somatorio += percentual
  print("Posto: ", i+1, " (", round(percentual), "% do total) => (acumulado de ", round(somatorio), "%)", sep="")
  
  if somatorio > 70 and posto_ideal_nao_descoberto:
    posto_ideal_nao_descoberto = False
    print("Posto ideal:", i+1)

Este seria o resultado: