Criando um preditor de proteínas de membrana com Python
Neste tutorial mostrarei como podemos criar uma ferramenta simples de bioinformática para determinar a localização sub-celular de uma proteína :D. Para isso, vamos utilizar a linguagem de programação Python e algumas bibliotecas para machine learning (Pandas, Sklearn e XGBoost), bioinformática (BioPython) e desenvolvimento web (Flask).
O que é localização sub-celular?
As proteínas podem atuar em diferentes compartimentos sub-celulares, como membrana, citoplasma, núcleo, espaço inter-membranas (no caso de bactérias Gram-negativas), bem como ser transportadas para fora da célula (secreção). Dependendo das características fisico-químicas do meio no qual a proteína deverá desempenhar a sua atividade, a sua composição de aminoácidos será diferente. Proteínas de membrana, por exemplo, tendem a ter uma predominância maior de aminoácidos mais hidrofóbico, enquanto que proteínas citoplasmáticas, de modo a manter a sua solubilidade, possuem aminoácidos mais hidrofílicos em sua estutura.
No nosso caso, criaremos um modelo que classificação as proteínas em apenas duas categorias, “citoplasmática” e “membrana”, a partir de sua composição de aminoácidos. Apesar de simplista, esta classificação é suficiente para ilustrar alguns dos principais conceitos sobre aprendizagem de máquina no contexto da bioinformática, e posteriormente poderemos extender ela de modo a abranger mais localizações, bem como incorporar dados adicionais para incrementar a sua confiabilidade :)
Instalando dependências
Neste projeto utilizaremos a linguagem Python 3 e algumas bibliotecas externas (também chamadas “pacotes”), sendo elas:
- XGBoost: pacote com algoritmos de regressão e classificação baseados em árvore de decisão com gradient boosting.
- BioPython: pacote para manipulação de dados biológicos.
- Sklearn: pacote para aprendizagem de máquina e manipulação de dados.
- Pandas: pacote para manipulação de dados tabulares.
- Flask: microframework para criação de aplicações web.
Para instalar elas é necessário ter o programa pip3
, o gerenciador de pacotes da linguagem Python 3. Caso você não tenha, pode instalar (no Ubuntu) com o comando:
$ sudo apt install python3-pip
Caso esteja utilizando o Windows, o pip3
provavelmente já estará instalado se estiver com a ultima versão do Python 3. Caso não esteja, atualize sua versão e durante a instalação marque para adicionar o Python e o Pip nas suas variáveis de ambiente.
Depois de instalar o pip3
, digitar o seguinte comando:
$ sudo pip3 install flask biopython sklearn xgboost pandas
Estrutura do projeto
Nosso projeto deverá conter a seguinte estrutura:
project/
scripts/
download.sh
preprocess.py
train.py
server.py
templates/
index.html
data/
Makefile
Os scripts Python que vamos utilizar ficarão na pastascripts
, sendo cada um responsável por uma etapa de análise. A etapa de preprocessamento dos dados (converter dados de sequências para um mais “amigável” para o algoritmo de predição) será realizada pelo scriptpreprocess.py
, enquanto que o treinamento do modelo preditivo será realizado pelo train.py
e a aplicação web será implementada no script server.py
.
Por padrão eu costumo organizar as etapas na ordem que elas devem ser executdas utilizando arquivos do tipo Makefile
. Neste caso, o conteúdo do arquivo fica assim:
preprocess:
python3 scripts/preprocess.pytrain:
python3 scripts/train.pyserver:
python3 scripts/server.py
Neste arquivo temos 3 etapas (regras) sendo definidas: preprocess
, train
e server
. Em cada uma especificada um script específico será executado.
Obtendo o dataset
Neste tutorial vamos precisar de dois arquivos no formato FASTA, um contento sequências de proteínas citoplasmáticas e outro contendo sequências de proteínas de membrana.
- Proteínas citoplasmáticas: baixe as sequências de proteínas citoplasmáticas disponíveis neste link e salve na pasta
data
com o nomecyt_proteins.fasta
. - Proteínas de membrana: baixe as sequências de proteínas de membrana disponíveis neste link e salve na pasta
data
com o nomemem_proteins.fasta
Preprocessando os dados
Nesta etapa vamos converter os dados de sequência para um formato tabular que poderá ser usado pelo algoritmo de predição. Esse processo de conversão é necessário pois para muitos algoritmos de aprendizagem de máquina, sequências de proteínas nada mais são que texto puro, enquanto que seu funcionamento, geralmente, é baseado no uso de inputs numéricos. Além disso, cada posição em uma sequência de proteínas reflete uma dimensão diferente, e como temos proteínas com tamanhos diferentes, os dados de nosso dataset possuem dimensionalidade variável, o seja, são dados não estruturados. Podemos resolver este problema de diversas formas, e uma das abordagem mais utilizadas para o problema de predição de localização sub-celular é a composição de aminoácidos.
A composição de aminoácidos é interessante por dois motivos: (1) reflete indiretamente características fisico-químicas das proteínas e (2) permite que todas sequências sejam representadas pelo mesmo número restrito de atributos (20, um para cada aminoácido). No script preprocess.py
, esta etapa é realizada da seguinte forma:
from Bio.SeqUtils.ProtParam import ProteinAnalysis
from Bio import SeqIO
import pandas as pddef get_protein_composition(sequence):
protein_analysis = ProteinAnalysis(sequence)
return protein_analysis.get_amino_acids_percent()cyt_proteins = SeqIO.parse(open("data/cyt_proteins.fasta"), 'fasta')
mem_proteins = SeqIO.parse(open("data/mem_proteins.fasta"), 'fasta')dataset_raw = []for protein in cyt_proteins:
protein_data = get_protein_composition(str(protein.seq))
protein_data['mem'] = False
dataset_raw.append(protein_data)for protein in mem_proteins:
protein_data = get_protein_composition(str(protein.seq))
protein_data['mem'] = True
dataset_raw.append(protein_data)pd.DataFrame(dataset_raw).to_csv("data/dataset.csv", index=False)
Todas análise da composição de aminoácidos é realizada dentro da função get_protein_composition
pelo objeto protein_analysis
que pertence à classe ProteinAnalysis
. Esta função recebe uma sequência e retorna um dicionário contendo a ocorrência de cada aminoácido em uma determinada proteína e é executada para ambos os datasets. Além da composição, adicionamos também no nosso novo dataset uma label para indicar a localização da proteína em questão (True
para proteínas de membrana e False
para proteínas citoplasmáticas). O novo dataset é então salvo na pasta data
com o nome dataset.csv
.
Treinando o modelo
O código do script que realiza o treinamento do modelo está apresentado abaixo.
from xgboost import XGBClassifier
from sklearn.model_selection import cross_val_score
import pandas as pd
import numpy as np
import pickle
import warningsdf_train = pd.read_csv("data/dataset.csv").sample(10000)
X = df_train.drop(['mem'], axis=1)
y = df_train['mem']xgb = XGBClassifier()
xgb.fit(X, y)
xgb_scores = cross_val_score(xgb, X, y, cv=10, scoring='roc_auc')
xgb_score = np.mean(xgb_scores)pickle.dump(xgb, open('data/model.pickle', 'wb'))print("XGB - AUC (ROC): ", xgb_score)
Neste caso, ao carregar os dados gerados na etapapreprocess
optei por selecionar apenas 10,000 entradas, e não o dataset completo. Esta escolha foi unicamente para tornar o processo de treinamento mais rápido, mas não tem maiores efeitos na acurácia do modelo final. Após carregar os dados para a variável df_train
, os atributos da sequência são passados para um variável chamada X
, e a label específica de cada sequência (True para proteínas de membrana, False para proteínas citoplasmástica) é passada para uma segunda variável chamada y
. Caso você já tenha lido sobre data science com Python deve estar reconhecimento este formato, pois é o mesmo usado em muitos materiais sobre Pandas e Sklearn.
Após carregar os dados para as variáveis, criamos um objeto chamado xgb
a partir da classe XGBClassifier
. Este objeto será o nosso modelo preditivo, e para ele são passadas as variáveis X e y através do método .fit(X,y)
. Com em outros algoritmos de aprendizagem de máquina com Python, o método fit()
da classe XGBClassifier
serve para iniciar o processo de treinamento.
Após o treinamento, utilizamos a função cross_val_score
para avaliar o modelo. Esta função realizar o processo de validação cruzada (cross-validade), onde o dataset é separado em n sub-sets, sendo n-1 usados para treinamento e 1 usado para a validação. Este processo é repetido n vezes, e para cada etapa é realizada a avaliação do modelo a partir de uma função de nosso interesse. Neste caso, escolhi a função área sob a curva ROC ( “roc_auc”) . Como para cada etapa é computador um valor de pontuação diferente, usei a função np.mean
para gerar uma média de todos os valores calculados.
Após avaliar o modelo, utilizamos as funções do pacote pickle
para salvar o modelo treinamento com o nome model.pickle
de modo que possamos utilizar ele futuramente.
Disponibilizando o modelo através de um webApp
Por fim, vamos criar uma pequena aplicação web que receba as sequências que queremos analisar. Para isso vamos utilizar o Flask, um microframework bem versátil do Python para desenvolvimento web. Não vou entrar em maiores detalhes sobre a estrutura do app pois isso foge do escopo deste tutorial, mas caso queiram estudar mais sobre esta ferramenta, seguem alguns links bem úteis :D [1, 2 e 3].
No contexto deste tutorial, as partes mais interessantes deste código são:
- Carregamento do modelo gerado no passo anterior com a função
pickle.load(...)
; - Conversão da sequência recebida para o mesmo formato usado no treinamento, sendo também usada a função
get_protein_composition
. - Uso do método
model.predict(..)
para predizer se a proteína é ou não de membrana (True
ouFalse
). - Uso do método
model.predict_proba(...)
para predizer a probabilidade da proteína ser de membrana (valor entre 0 e 1).
from flask import Flask, request, redirect, render_template
from Bio.SeqUtils.ProtParam import ProteinAnalysis
import pickle
import pandas as pdapp = Flask(__name__)
model = pickle.load(open("data/model.pickle", "rb"))def get_protein_composition(sequence):
protein_analysis = ProteinAnalysis(sequence)
return protein_analysis.get_amino_acids_percent()@app.route("/")
def index():
return render_template('index.html')@app.route("/predict", methods=['GET', 'POST'])
def predict():
if request.method == "GET":
return redirect("/")
sequence = request.form['sequence']
protein_data = get_protein_composition(sequence)
dataframe = pd.DataFrame([protein_data])
prediction = model.predict(dataframe)[0]
probability = model.predict_proba(dataframe)[0][1]if prediction:
location = "Membrane"
else:
location = "Cytoplasm"return render_template('index.html',
result={"location": location,
"probability": probability,
"sequence": sequence})if __name__ == "__main__":
app.run(host="0.0.0.0", debug=True)
Para que este código funcione, precisamos criar um arquivo HTML que será usado como template pelo Flask (mais especificamente, pelo Jinja).
Nosso arquivo template deverá ter o nome index.html
e deve estar localizado em uma pasta com o nome templates
dentro da mesma pasta onde está o script da nossa aplicação web (server.py
).
<!DOCTYPE html>
<html lang="en">
<body>
<h1>Membrane Protein Predictor</h1>
<form action="/predict" method="POST">
<textarea name="sequence" placeholder=">sequence ..."
style="width: 400px; height: 100px"></textarea><br>
<input type="submit" value="submit">
</form>
<hr>
{% if result %}
<h2>Result of Analysis:</h2>
<pre>{{result['sequence']}}</pre>
Protein Location: {{result['location']}} <br>
Probability (membrane): {{result['probability']}} <br>
{% endif %}
</body>
</html>
Executando tudo
Agora, basta executar os scripts com o Makefile
! :D
Preprocessamento
Comando: make preprocess
$ make preprocess
$ ls data/
cyt_proteins.fasta dataset.csv mem_proteins.fasta
Treinamento
Comando: make train
$ make train
XGB - AUC (ROC): 0.9131659835214482
Executando a aplicação
Comando: make server
$ make server
* Serving Flask app "server" (lazy loading)
* Environment: production
WARNING: This is a development server. Do not use it in a production deployment.
Use a production WSGI server instead.
* Debug mode: on
* Running on http://0.0.0.0:5000/ (Press CTRL+C to quit)
* Restarting with stat
* Debugger is active!
* Debugger PIN: 252-910-284
Agora só acessar o endereço localhost:5000
no navegador para conferir a ferramenta rodando! :D
Agora, basta colocar uma sequência de proteína para ser analizar ali na text area e clicar em submit!
Para testar, podemos utilizar uma sequência conhecida de proteína de membrana, como a bomba de efluxo de tolueno de Pseudomonas putida (UniProt: Q93PU3) (Obs: copie a sequência sem o header do FASTA).
Feito! Criamos nosso primeiro preditor de localização sub-celular! Ainda podemos melhorar a sua confiabilidade incorporando mais informações, bem como aumentar o número de localizações possíveis, mas de modo geral a estrutura será esta mesma!