Identificador de Hidrocarbonetos Alifáticos Saturados pelo padrão InChI¶
Visão Computacional e Processamento de Imagens¶
Jhonata Santana Louzada de Aguiar - 11201811899¶
Prof. Francisco Zampirolli - 2019q3¶
Introdução¶
Motivação¶
Atualmente estou cursando a disciplina de Bioquímica e a professora tem certa facilidade com a identificação de compostos, quem não possui interesse na área de química passa certa dificuldade para resolver problemas do tipo, ou até mesmo ao tentar resolver alguns exercícios envolvendo a estrutura molecular de um composto, pois não possuem esta facilidade. Convenhamos que para quem não possui interesse na área acaba sendo um desperdício de memória, pois querendo ou não é uma baita carga de regras a serem decoradas. Há diversas ferramentas que auxiliam na resolução de problemas matemáticos, por que não há uma conhecida e eficiente para solução de problemas químicos? Eis que após resolver um problema matemático no Photomath surge a ideia de tentar criar uma ferramenta semelhante, mas que inicialmente sua utilidade seria apenas a de identificar o nome do composto representado em uma imagem, afinal, possuindo esta informação a resolução de problemas pode ser facilitada.
Projeto¶
Análise de imagens de moléculas de hidrocarbonetos para gerar o código no padrão InChI representante desta para a consulta de informações no banco de dados da ChemSpider, ou até mesmo a disponibilização deste código para a consulta no Google e outras ferramentas.
Ferramentas¶
- Python
- Tesseract OCR
- ChemSpider
- PubChem(apenas para teste e verificações)
O padrão InChI da IUPAC¶
Sinceramente não possuia nenhuma referência a respeito deste padrão e muitas pessoas, inclusive professores da área também desconhecem tal método. "O Identificador Químico Internacional da IUPAC (InChI - International Chemical Identifier) é um identificador textual para substâncias químicas, com o objetivo de estabelecer uma maneira padrão de descrever informações de moléculas e facilitar sua pesquisa." - Wikipédia. Este padrão ganhou mais força, embora quase ninguém o conheça, em 2011 e sua proposta é simplesmente poder descrever a estrutura de uma molécula em apenas uma linha, esta por sua vez é dividida em "camadas"(/) onde cada uma é responsável por armazenar certas informações, solucionando problemas presentes em outros identificadores. Para este projeto trabalharemos apenas com 4 camadas:
- a responsável pela versão do identificador;
- a camada responsável pela fórmula geral da molécula;
- a camada responsável por informar a estrutura/ligações da molécula;
- e por fim a camada responsável pela disposição dos hidrogênios.
O propano(C3H8) possui o seguinte código: InChI=1S/C3H8/c1-3-2/h3H2,1-2H3
plt.figure(figsize=(15, 15))
plt.axis('off')
plt.imshow(cv.cvtColor(cv.imread("propano.jpg"), cv.COLOR_RGB2BGR))
plt.show()
ChemSpider¶
O ChemSpider é um banco de dados que conta com mais de 67 milhões de estruturas de moléculas quem podem ser pesquisadas desde sua nomenclatura padrão literal até códigos mais complexos como o InChI. Há outro banco de dados o do PubChem, por exemplo, que pemite realizar as mesmas pesquisas, porém apenas o ChemSpider, pelo que encontrei até o momento, possui uma API que permite acesso a este banco de forma gratuita, embora sejam permitidas apenas 1000 consultas por mês, o que é mais do que o suficiente para elaboração do projeto.
Por que o InCHI?¶
A escolha do padrão Inchi foi feita justamente pela proposta do padrão, ao qual até o momento é o que não possui lá muitas falhas encontradas, sem contar com o suporte da Google para tal, onde podemos realizar uma pesquisa através do mecanismo de busca do Google utilizando este codificador e obteremos como primeiros resultados, inclusive as imagens, conteúdos pertinentes a este composto.
plt.figure(figsize=(15, 15))
plt.axis('off')
plt.imshow(cv.cvtColor(cv.imread("googles.png"), cv.COLOR_RGB2BGR))
plt.figure(figsize=(15, 15))
plt.axis('off')
plt.imshow(cv.cvtColor(cv.imread("googlei.png"), cv.COLOR_RGB2BGR))
plt.show()
Mão na massa¶
Começamos com as importações
import pytesseract
import cv2 as cv
import numpy as np
import urllib.request
import matplotlib.pyplot as plt
from PIL import Image
from operator import itemgetter
from chemspipy import ChemSpider
# Redefinição do PATH onde se encontra o Tesseract (variável para cada máquina)
pytesseract.pytesseract.tesseract_cmd = r'C:\Program Files (x86)\Tesseract-OCR\tesseract.exe'
Em seguida definimos algumas funções
# Função responsável por adicionar bordas em imagens com apenas 1 canal
def addBordaGray(image, borda, val):
img = (np.ones((image.shape[0] + borda*2, image.shape[1] + borda*2), dtype=int))*val
img[borda:image.shape[0]+borda, borda:image.shape[1]+borda] = image.copy()
return img
# Função responsável por atribuir o próximo "nó" na construção do código InChI
def nextnode(node, pos):
connection = atoms[int(node["lig"][pos][1:]) - 1]
nodes[connection["peso"] - 1]["lig"].remove(node["elemento"] + str(node["id"]))
node = nodes[connection["peso"] - 1]
return node
# Função responsável pela construção da camada de estrutura do código InChI
def connect(node, inchi):
if len(node["lig"]) != 0:
inchi += str(node["peso"])
if len(node["lig"]) == 1:
inchi += '-'
inchi = connect(nextnode(node, 0), inchi)
elif len(node["lig"]) == 2:
inchi += "("
inchi = connect(nextnode(node, 0), inchi)
inchi += ")"
inchi = connect(nextnode(node, 1), inchi)
elif len(node["lig"]) == 3:
inchi += "("
inchi = connect(nextnode(node, 0), inchi)
inchi += ","
inchi = connect(nextnode(node, 1), inchi)
inchi += ")"
inchi = connect(nextnode(node, 2), inchi)
else:
inchi += str(node["peso"])
return inchi
# Função adicional exibição dos "pesos" de cada "nó" para verificação da
# execução dos procedimentos, assim como o algoritmo de Morgan
def mostraPeso(nodes, atoms, image):
img = image.copy()
for box in nodes:
cv.putText(img, str(box["peso"]), (atoms[box["id"] - 1]["centro"][0]-10, atoms[box["id"] - 1]["centro"][1]), cv.FONT_HERSHEY_SIMPLEX, 1, 0, 2,
cv.LINE_AA)
plt.figure(figsize=(10, 10))
plt.axis('off')
plt.imshow(img, "gray")
plt.show()
Recebemos a imagem que será analisada e a preparamos para as futuras análises
image = cv.imread("eter.jpg")
image = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
# Essas duas imagens a seguir serão utilizadas futuramente para análise, pois a variável
# original será utilizada em outros processos irreversíveis
borda = 10
imgaux = addBordaGray(image, borda, 255)
connection = imgaux.copy()
# image = cv.GaussianBlur(image, (3,3), 0)
dump, image = cv.threshold(image, 10, 255, cv.THRESH_BINARY_INV+cv.THRESH_OTSU)
contours, hierarquia = cv.findContours(image, cv.RETR_TREE, cv.CHAIN_APPROX_SIMPLE)
for area in contours:
cv.drawContours(image, [area], 0, 200, -1)
Em seguida procuramos os atómos e ligações presentes na mólecula da imagem por meio de seus contornos, então obtemos sua bounding box
contours, hierarquia = cv.findContours(image, cv.RETR_TREE, cv.CHAIN_APPROX_SIMPLE)
boxes = []
for obj in contours:
minx, maxx, miny, maxy = len(image[0]), 0, len(image), 0
for point in obj:
if point[0, 0] < minx:
minx = point[0, 0]
elif point[0, 0] > maxx:
maxx = point[0, 0]
if point[0, 1] < miny:
miny = point[0, 1]
elif point[0, 1] > maxy:
maxy = point[0, 1]
boxes.append((minx, maxx, miny, maxy))
Agora que possuímos as bounding boxes de cada átomo e ligação devemos diferenciá-los
ligacoes, atoms = [], []
for box in boxes:
height = box[3] - box[2]
width = box[1] - box[0]
cx = int(width / 2)
cy = int(height / 2)
centro = (box[0] + cx, box[2] + cy)
if 0 <= width < height / 2 or 0 <= height < width / 2:
if 0 <= width < height / 2:
tipo = 'v'
elif 0 <= height < width / 2:
tipo = 'h'
ligacao = {"id": len(ligacoes) + 1,
"tipo": tipo,
"xi": box[0],
"xf": box[1],
"yi": box[2],
"yf": box[3],
"width": width,
"height": height,
"centro": centro,
"ligacao1": 0,
"ligacao2": 0}
ligacoes.append(ligacao)
else:
letter = {"id": len(atoms) + 1,
"letra": 0,
"ligacao": [],
"LH": 0,
"LC": 0,
"peso": 1,
"xi": box[0],
"xf": box[1],
"yi": box[2],
"yf": box[3],
"width": width,
"height": height,
"centro": centro,
"d1":0}
atoms.append(letter)
Agora que sabemos quais são os átomos e onde estão, vamos identificar cada um deles formando uma palavra com cada átomo
# utilizamos o átomo com caracter mais alto como para a altura da palavra.
maxheight = max(atom["height"] for atom in atoms)
margin = 60
word = (np.ones((maxheight + margin, margin), dtype=int))*255
erro = 3
for box in atoms:
letra = (np.ones((maxheight + margin, box["width"] + 8), dtype=int))*255
# Recortamos a letra da imagem
letra[int(margin/2) - erro:box["height"] + int(margin/2) + erro, 4 - erro:box["width"] + 4 + erro] = imgaux[box["yi"] + borda - erro:box["yf"] + borda + erro, box["xi"] + borda - erro:box["xf"] + borda + erro].copy()
# Concatenamos na palavra
word = np.concatenate((word, letra), axis=1)
end = (np.ones((maxheight + margin, margin), dtype=int))*255
word = np.concatenate((word, end), axis=1)
word = word.astype("uint8")
word = cv.GaussianBlur(word, (3,3), 0)
dump, word = cv.threshold(word, 50, 255, cv.THRESH_BINARY+cv.THRESH_OTSU)
plt.figure(figsize=(15, 15))
plt.axis('off')
plt.imshow(word, "gray")
plt.show()
#Analisa o que está escrito e retorna os elementos em uma string
word = Image.fromarray(word)
txt = pytesseract.image_to_string(word)
txt = txt.upper()
txt = txt.replace('0', 'O')
txt = txt.replace('4', 'H')
txt = txt.replace(' ', '')
txt = txt.replace('.', '')
txt = txt.replace('-', '')
txt = txt.replace(',', '')
txt = txt.replace(';', '')
txt = txt.replace('=', '')
txt = txt.replace('_', '')
txt = txt.replace('"', '')
txt = txt.replace("''", '')
txt = txt.replace('´', '')
txt = txt.replace('`', '')
Eis que obtemos a seguinte string:
print(txt)
Cada letra representa o átomo de cada bounding box, então a seguir iremos atribuir ao dicionário, assim como já criar um dicionário para armazenar os átomos existentes e suas respectivas quantidades
elements = {}
for i in range(len(txt)):
atoms[i]["letra"] = txt[i]
if txt[i] not in elements:
elements[txt[i]] = 0
if txt[i] in elements:
elements[txt[i]] +=1
Em seguida já é possível obter a fórmula geral da molécula, responsável pela 2ª camada do código InChI
formula = ""
#Inicia a fórmula da molécula conforme as regras, C, H e o restante dos elementos em ordem alfabética.
if 'C' in elements:
formula += 'C'
if elements['C'] > 1:
formula += str(elements['C'])
if 'H' in elements:
formula += 'H'
if elements['H'] > 1:
formula += str(elements['H'])
#Ordena os elementos em ordem alfabética.
elements = sorted(elements.items(), key=itemgetter(1))
#Finaliza o registro da molécula com os elementos diferentes de C e H.
for i in elements:
if i[0] != 'C' and i[0] != 'H':
formula += str(i[0])
if i[1] > 1:
formula += str(i[1])
Agora iremos obter as informações necessárias para a 3ª camada
edges = []
#Identifica os participantes de cada ligação.
for lig in ligacoes:
lig1 = 0
lig2 = 0
if lig["tipo"] == 'v':
topdif = 1000
bottomdif = 1000
for ele in atoms:
if ele["xi"] < lig["centro"][0] < ele["xf"]:
if ele["yf"] < lig["yi"] and 0 < (lig["yi"] - ele["yf"]) < topdif:
topdif = lig["yi"] - ele["yf"]
lig1 = ele["letra"] + str(ele["id"])
elif ele["yi"] > lig["yf"] and 0 < (ele["yi"] - lig["yf"]) < bottomdif:
bottomdif = ele["yi"] - lig["yf"]
lig2 = ele["letra"] + str(ele["id"])
else:
leftdif = 1000
rightdif = 1000
for ele in atoms:
if ele["yi"] < lig["centro"][1] < ele["yf"]:
if ele["xf"] < lig["xi"] and 0 < (lig["xi"] - ele["xf"]) < leftdif:
leftdif = lig["xi"] - ele["xf"]
lig1 = ele["letra"] + str(ele["id"])
elif ele["xi"] > lig["xf"] and 0 < (ele["xi"] - lig["xf"]) < rightdif:
rightdif = ele["xi"] - lig["xf"]
lig2 = ele["letra"] + str(ele["id"])
# Interpreta as ligações e as atribui para cada elemento com seu respectivo id.
lig["ligacao1"] = lig1
lig["ligacao2"] = lig2
ponta1 = int(lig1[1:]) - 1
ponta2 = int(lig2[1:]) - 1
if lig2[0] != 'H':
atoms[ponta1]["ligacao"].append(lig2)
if lig1[0] != 'H':
atoms[ponta2]["ligacao"].append(lig1)
# Atribui a quantidade de ligações com Hidrogênios e sem.
if lig2[0] == 'H':
atoms[ponta1]["LH"] += 1
else:
atoms[ponta1]["LC"] += 1
if lig1[0] == 'H':
atoms[ponta2]["LH"] += 1
else:
atoms[ponta2]["LC"] += 1
edges.append((lig1, lig2))
Agora que já sabemos quais átomos interagem entre si vamos refinar estas informações e separa apenas o que interessa para realizar a próximas operações, que são os átomos diferentes de Hidrogênio, ao qual chamaremos de nós
nodes = []
for i in edges:
if i[0][0] != 'H' and i[1][0] != 'H':
for j in i:
if j not in nodes:
nodes.append(j)
# Obtemos as informações destes átomos e eliminamos as ligações de Hidrogênio da lista de ligações
if len(nodes) != 0:
for i in range(len(nodes)):
nodes[i] = {"id":int(nodes[i][1:]),
"elemento":nodes[i][0],
"peso":1,
"lig":sorted(atoms[int(nodes[i][1:]) - 1]["ligacao"]),
"lh":atoms[int(nodes[i][1:]) - 1]["LH"],
"leftd":atoms[int(nodes[i][1:])]["xi"],
"d1":0,
"cc":0}
Agora executamos o algoritmo de Morgan, essencial para esta 3ª camada, ao qual irá nos orientar na estrutura da molécula
if len(nodes) != 0:
x = 1
v = 0
aux = 0
mostraPeso(nodes, atoms, connection)
while x != 0:
qtdv = []
for j in range(len(nodes)):
nodes[j]["peso"] = 0
for k in nodes[j]["lig"]:
nodes[j]["peso"] += atoms[int(k[1:]) - 1]["peso"]
for m in nodes:
atoms[m["id"] - 1]["peso"] = m["peso"]
if m["peso"] not in qtdv:
qtdv.append(m["peso"])
if len(qtdv) > v:
v = len(qtdv)
aux = 0
elif len(qtdv) == v and aux != 2:
aux += 1
if len(qtdv) == v and aux == 2:
x = -1
x += 1
mostraPeso(nodes, atoms, connection)
Agora ordenamos os pesos a partir de 1
if len(nodes) != 0:
# Ordena os nós de acordo com seus pesos e classifica o primeiro como 1 para dar início
nodes = sorted(nodes, key= itemgetter("elemento", "peso", "leftd"))
nodes[0]["peso"] = 1
atoms[nodes[0]["id"] - 1]["peso"] = 1
# Atribui a distãncia do primeiro nó, utiizado como parâmetro de desempate
for i in range(1, len(nodes)):
nodes[i]["d1"] = (((atoms[nodes[i]["id"] - 1]["centro"][0] - atoms[nodes[0]["id"]]["centro"][0])**2) + ((atoms[nodes[i]["id"] - 1]["centro"][1] - atoms[nodes[0]["id"]]["centro"][1])**2))**(1/2)
# Ordena os nós de acordo com seus pesos e distância do átomo inicial
nodes = sorted(nodes, key= itemgetter("peso", "d1"))
# Atribui os pesos reais para obter o código InChI
for i in range(1, len(nodes)):
nodes[i]["peso"] = i + 1
atoms[nodes[i]["id"] - 1]["peso"] = nodes[i]["peso"]
Agora que sabemos o peso de cada átomo devemos ordenar a lista de ligações dos nós de acordo com o peso de cada um deles
if len(nodes) != 0:
for i in range(len(nodes)):
vet = []
ligs = nodes[i]["lig"].copy()
for j in range(len(nodes[i]["lig"])):
vet.append(atoms[int(nodes[i]["lig"][j][1:]) - 1]["peso"])
compare = sorted(vet.copy())
nodes[i]["lig"] = []
for k in compare:
pos = vet.index(k)
nodes[i]["lig"].append(ligs[pos])
Agora vem a cereja do bolo, gerar a 3ª camada
inchi = "InChI=1S/" # A primeira camada é simples, 1S, versão do InChI
inchi += formula # Incrementamos a 2ª camada
# Geramos a 3ª camada
if len(nodes) != 0:
inchi += "/c"
node = nodes[0]
inchi = connect(node, inchi)
print(inchi)
Como estamos tratando de hidrocarbonetos saturados, claramente há hidrogênios, então devemos formar a 4ª camada
# Geramos a 4ª camada
if len(nodes) != 0:
inchi += "/h"
x = 0
while elements[x][0] != 'H':
x += 1
posh = [0] * (elements[x][1] + 1)
for i in range(len(posh)):
posh[i] = [0]
for i in nodes:
posh[i["lh"]][0] += 1
posh[i["lh"]].append(i["peso"])
aux = 0
for i in range(1, len(posh)):
if posh[i][0] != 0:
if aux == 1:
inchi += ','
if len(posh[i]) == 2:
inchi += str(posh[i][1]) + 'H'
aux = 1
elif posh[i][0] == (max(posh[i][1:]) - min(posh[i][1:]) + 1):
inchi += str(min(posh[i][1:])) + '-' + str(max(posh[i][1:])) + 'H'
aux = 1
if i != 1:
inchi += str(i)
print(inchi)
Agora que conseguimos gerar o código InChI, basta utilizar a bibloteca chemspipy para pesquisar informações sobre esta molécula
cs = ChemSpider('2bEGovZUrX2htIddZBOZFv4DUpzDvwqo') # Utilizo a key da API que obtive
queryID = cs.filter_inchi(inchi)
queryResults = cs.filter_results(queryID)
if len(queryResults) != 0:
molecula = cs.get_compound(queryResults[0])
print(f"""
Nome Comum: {molecula.common_name}
Fórmula Molecular: {formula}
InChI: {inchi}
InChIKey: {molecula.inchikey}
SMILES: {molecula.smiles}
Peso Molecular: {molecula.molecular_weight}
Massa Média: {molecula.average_mass}
Massa Monoisotópica: {molecula.monoisotopic_mass}
Massa Nominal: {molecula.nominal_mass}
""")
imgmol = Image.open(urllib.request.urlopen(molecula.image_url))
imgmol = cv.cvtColor(np.array(imgmol), cv.COLOR_RGB2BGR)
plt.figure(figsize = (5,5))
plt.imshow(imgmol)
plt.show()
else:
print("Não foi possível identificar esta molécula, tente outra imagem.")
Podemos encontrar diversos link atrelados a molécula identificada
if len(queryResults) != 0:
print("Referências Externas:")
links = []
for ref in molecula.external_references[:10]:
if ref["externalUrl"] not in links and ref["externalUrl"] != '':
links.append(ref["externalUrl"])
print(ref["externalUrl"])
print("...")
O programa está apto a receber outros elementos, desde que estejam saturados e devem possuir suas ligações explícitas. No entanto como veremos a seguir, possuímos algumas limitações decorrentes da ineficiência do algoritmo de Morgan e a ausência de regras para desempate para o mesmo.
molIdentifier("isopentano.png")
Como podemos observar ele não identifica a molécula do isopentano, mas caso seja feita a alteração de um dos parâmetros da codificação responsável pelo algoritmo de Morgan, permitindo mais uma execução do método, este composto será identificado perfeitamente, faremos um teste a seguir alterando este parâmetro
molIdentifier("isopentano.png", 3)
Como vimos, uma repetição de diferença pode alterar o código InChI, essa diferença ocorre pois há certos empates no algoritmo de Morgan que exige um conhecimento e a implementação de métodos específicos para o desempate, sendo assim, acredito que seja fundamental a implementação de uma inteligência artificial para auxílio na tomada de decisões e orientação a objeto para melhorar a performance e organização do algoritmo