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.

Transportadores de glicose são um exemplo de proteínas de membrana (fonte da imagem: PDB).

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.py
train:
    python3 scripts/train.py
server:
    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 nome cyt_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 nome mem_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 pd
def 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 warnings
df_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 ou False).
• 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 pd
app = 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!