O objectivo desta fase era, dadas as proteínas que retirámos do ficheiro genbank, procurar sequências homólogas com a ferramenta BLAST. Esta fase foi iterada várias vezes à procura da melhor estratégia, e no final acabámos com várias alternativas para a execução do BLAST, que podem ser encontradas no módulo util.blast. As alternativas são:
Os resultados desta fase podem ser encontrados aqui.
Esta função assume que o BLAST está instalado localmente, e recebe como parâmetro a base de dados que deve usar. O primeiro argumento desta função, tag_to_files, é um dicionário no qual as chaves são os locus_tag de cada gene, e os valores são pares em que a primeira componente é o ficheiro, onde se encontra a sequência de aminoácidos para ser usada na query, (in_file), e a segunda componente é o ficheiro onde deve ser guardado o resultado do BLAST (out_file).
def local_blastp(tag_to_files, db):
"""
Corre o blast localmente.
"""
for tag in tag_to_files:
(in_file, out_file) = tag_to_files[tag]
blastp_cline = NcbiblastpCommandline(
query=in_file,
db=db,
evalue=10,
outfmt=5,
out=out_file
)
blastp_cline()
Para evitar a instalação do BLAST, o download de bases de dados e a sua adição como base de dados a ser usada pelo BLAST, criámos uma imagem Docker (vitorenesduarte/swissprot_blast) com a base de dados UniProtKB - Swiss-Prot retirada do site da UniProt.
def docker_blastp(directory, db):
"""
Corre o blast numa instância docker.
"""
cmd = "docker run -e QUERY_DIR=" + directory \
+ " -e DB=" + db \
+ " -v $PWD/" + directory + ":/" + directory \
+ " -ti vitorenesduarte/swissprot_blast"
p = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
p.wait()
A terceira opção não corre o BLAST localmente, mas utiliza o serviço disponibilizado pelo site ExPASy.
def expasy_blastp(tag_to_files):
"""
Corre o blast na expasy.
"""
for tag in tag_to_files:
(in_file, out_file) = tag_to_files[tag]
query = rw.read_file(in_file)
blast_result = www.expasy_blast(query)
rw.write_file(blast_result, out_file)
Esta função usa uma outra função definida no módulo util.www que corre o BLAST contra a base de dados UniProtKB = Swiss-Prot + TrEMBL, usando o serviço da ExPASy.
def expasy_blast(query):
"""
Corre o blast na expasy.
"""
url = "http://web.expasy.org/cgi-bin/blast/blast.pl"
data = {
"format": "xml",
"prot_db1": "UniProtKB",
"matrix": "BLOSUM62",
"showsc": 50, # best scoring sequences
"showal": 50, # best alignments
"seq": query
}
response = requests.post(
url,
data=data
)
return response.text
Esta última opção permitiu-nos encontrar sequências com maior homologia do que as que estávamos a encontrar, utilizando apenas a Swiss-Prot localmente.
A função que engloba estas três opções recebe como argumento uma lista de pares em que a primeira componente é a locus_tag da proteína, e a segunda componente é a respectiva sequência de aminoácidos. Recebe ainda qual o tipo de BLAST, local, docker ou expasy, a utilizar.
def blastp(tags_and_proteins, type="local"):
"""
Corre o blast para a proteínas passadas como argumento.
O argumento type é opcional e pode ter dois valores:
- local
- docker
- expasy
É retornado um dicionário com os resultados.
"""
directory = ".query_dir"
tag_to_files = write_queries_to_dir(tags_and_proteins, directory)
# correr o blast
if type == "local":
local_blastp(tag_to_files, "swissprot")
elif type == "docker":
docker_blastp(directory, "swissprot")
elif type == "expasy":
expasy_blastp(tag_to_files)
else:
raise Exception("Unsupported type: " + type)
blast_results = extract_blast_info(tag_to_files, type)
return blast_results
Esta função usa duas outras funções:
write_queries_to_dir que recebe uma lista de pares, locus_tag e sequência de aminoácidos protein, e retorna um dicionário que mapeia locus_tag para pares (in_file, out_file).def write_queries_to_dir(tags_and_proteins, directory):
"""
Grava as proteínas passadas como argumento na
diretoria também passada como argumento.
"""
tag_to_files = {}
# Apagar a diretoria e criar uma nova
if os.path.exists(directory):
shutil.rmtree(directory)
os.makedirs(directory)
for (tag, protein) in tags_and_proteins:
# Gravar a proteína num ficheiro
in_file = directory + "/" + fasta_it(tag)
rw.write_file(protein, in_file)
# registar esta informação num dicionário
out_file = directory + "/" + xml_it(fasta_it(tag))
tag_to_files[tag] = (in_file, out_file)
return tag_to_files
extract_blast_info que dado este dicionário retorna um outro dicionário que mapeia locus_tag para a lista de resultados do BLAST. Cada um destes resultados é um dicionário com:def extract_blast_info(tag_to_files, type):
"""
Extrai a informação que necessitamos dos resultados do blast.
- uniprot_id
- evalue
- score
- identity
"""
blast_results = {}
for tag in tag_to_files:
(_, out_file) = tag_to_files[tag]
### Alguns resultados do blast, ao serem lidos de disco
### na função rw.read_blast
### estavam a dar erro ao fazer parsing nesta linha.
### Como não necessitamos deste resultado,
### optamos por removê-la.
regex = ".*Hsp_midline.*"
rw.filter_file(regex, out_file)
handle = rw.read_blast(out_file)
result = []
for a in handle.alignments:
# extrair o uniprot id
if type in ["local", "docker"]:
# nos blast locais, o uniprot_id está no hit_def
uniprot_id = a.hit_def.split("|")[1]
elif type == "expasy":
# nos blast expasy, o uniprot_id está no hit_id
uniprot_id = a.hit_id.split("|")[1]
# escolher sempre o primeiro hsp
hsp = a.hsps[0]
evalue = hsp.expect
score = hsp.score
identities = hsp.identities
align_length = hsp.align_length
identity = (identities * 100) / align_length
hit = {}
hit["uniprot_id"] = uniprot_id
hit["evalue"] = evalue
hit["score"] = score
hit["identity"] = identity
result.append(hit)
blast_results[tag] = result
return blast_results
Nesta fase ainda precisávamos de mais informação para além daquela que é possível extrair dos resultados do BLAST, que foi conseguida utilizando um serviço do site UniProt descrito aqui.
Esta informação foi adicionada a cada um dos resultados provenientes do BLAST.
def add_info_to_blast_results(blast_results, uniprots):
"""
Esta função adiciona aos blast results informação
retirada do site da uniprot.
"""
for tag in blast_results:
for i in range(len(blast_results[tag])):
uniprot_id = blast_results[tag][i]["uniprot_id"]
# esta proteína aparece nos resultados do blast como
# Q8TC84-2 mas a uniprot retorna a sua informação
# como Q8TC84
if uniprot_id == "Q8TC84-2":
uniprot_id = "Q8TC84"
properties = [
"status",
"organism",
"translation",
"molecular_functions"
]
for p in properties:
value = uniprots[uniprot_id][p]
blast_results[tag][i][p] = value
return blast_results
Por fim, automatizámos o processo de inferência da função, dados os resultados do BLAST. Uma função é considerada provável se aparecer em pelo menos 50% dos resultados do BLAST. Este inferência foi feita tanto para os 10 primeiros resultados como para os resultados todos.
def infer_function(blast_results):
"""
Infere qual as funções mais prováveis tendo em consideração
os primeiros 10 resultados do blast,
e quais as funções mais prováveis tendo em consideração
todos os resultados do blast.
"""
all = {}
for tag in sorted(blast_results.keys()):
results = blast_results[tag]
# inferir função para os primeiros 10 resultados
inferred_top_ten = infer(results[0:10])
# inferir função para todos os resultados
inferred_all = infer(results)
all[tag] = {}
all[tag]["top_ten"] = inferred_top_ten
all[tag]["all"] = inferred_all
return all
def infer(results):
"""
Esta função tenta inferir a função das proteínas baseando-se
numa lista de resultados do blast.
"""
inferred = []
leaderboard = defaultdict(int)
for result in results:
functions = result["molecular_functions"]
# atualizar a leaderboard
for function in result["molecular_functions"]:
leaderboard[function] += 1
# As funções que aparecem em pelo menos 50% dos resultados
# são consideradas potenciais funções
min = len(results) / 2
for function in leaderboard:
if leaderboard[function] >= min:
inferred.append(function)
return sorted(inferred)
Todos estes passos podem ser encontrados no ficheiro second.py e resumem-se a:
import util.rw as rw
import util.blast as blast
import util.www as www
import util.util as util
ncbi_json_path = ".ncbi.json"
blast_results_json_path = ".blast_results.json"
inferred_json_path = ".inferred.json"
dictionary = rw.read_json(ncbi_json_path)
tags_and_proteins = util.get_tags_and_proteins(dictionary)
blast_results = blast.blastp(
tags_and_proteins,
"expasy"
)
uniprot_ids = set()
for tag in blast_results:
for result in blast_results[tag]:
uniprot_ids.add(result["uniprot_id"])
uniprot_ids = list(uniprot_ids)
uniprots = www.fetch_uniprots(uniprot_ids)
blast_results = add_info_to_blast_results(
blast_results,
uniprots
)
rw.write_json(blast_results, blast_results_json_path)
inferred = infer_function(blast_results)
rw.write_json(inferred, inferred_json_path)