Os resultados desta fase podem ser encontrados aqui.
Começamos por aceder ao NCBI para fazer download de um ficheiro no formato genbank com informação correspondente à nossa zona do genoma, da posição 270001 a 505535.
O primeiro passo foi descobrir qual o genbank_id associado ao accession number, NC_002942.5, da estirpe em estudo.
Para tal, usamos a função Entrez.search.
Tendo o genbank_id, obtivemos o ficheiro genbank, com recurso à função Entrez.efetch.
def fetch_genbank(start, end):
"""
Procura no NCBI pelo Accession NC_002942.5 e
faz download de toda a informação em formato genbank,
retornando o respectivo record com as features
entre a posição start e end.
"""
Entrez.email = "vitorenesduarte@gmail.com"
handle = Entrez.esearch(db="nucleotide", term="NC_002942.5")
record = Entrez.read(handle)
genbank_id = record["IdList"][0]
handle.close()
handle = Entrez.efetch(db="nucleotide",
rettype="gbwithparts",
retmod="text",
id=genbank_id,
seq_start=start,
seq_stop=end)
record = SeqIO.read(handle, "gb")
handle.close()
return record
Esta função foi adicionada ao módulo util.www.
Deste record extraímos as features gene, CDS, tRNA e rRNA.
def extract_features(record):
"""
Das várias features encontradas no ficheiro genbank, extraímos:
- gene
- CDS
- tRNA
- rRNA
Esta função retorna uma lista apenas com as features
dos tipos indicados acima.
"""
result = []
features = record.features
feature_count = len(features)
types = ["gene", "CDS", "tRNA", "rRNA"]
for i in range(feature_count):
feature = features[i]
should_extract = feature.type in types
if should_extract:
result.append(feature)
return result
De seguida, convertemos esta lista de features num dicionário. A cada uma retirámos:
Além disso, também ficamos com informação sobre a localização de cada gene, cadeias em que se encontra, número de aminoácidos da proteína codificada e tipo de feature (tRNA, rRNA e mRNA).
def features_to_dictionary(start, features):
"""
Para cada uma das features, extraimos a seguintes propriedades:
- db_xref
- EC_number
- function
- gene
- gene_synonym
- locus_tag
- note
- product
- protein_id
- translation
Também extraímos a localização:
- start
- end
- strand
"""
dictionary = {}
properties = ["db_xref",
"EC_number",
"function",
"gene",
"gene_synonym",
"note",
"product",
"protein_id",
"translation"]
for feature in features:
tag = feature.qualifiers["locus_tag"][0]
if not tag in dictionary:
# se esta tag ainda não existe,
# criar um dicionário vazio para ela
dictionary[tag] = {}
location = feature.location
dictionary[tag]["location"] = {}
dictionary[tag]["location"]["start"] = location._start + start
dictionary[tag]["location"]["end"] = location._end + start - 1
dictionary[tag]["location"]["strand"] = location._strand
if feature.type in ["tRNA", "rRNA"]:
# Taggar as features com estes dois tipos
dictionary[tag]["type"] = feature.type
else:
# se não, taggar com "mRNA"
dictionary[tag]["type"] = "mRNA"
for prop in properties:
if prop in feature.qualifiers:
value = feature.qualifiers[prop][0]
if prop == "db_xref":
# renomear esta propriedade e remover "GeneID:"
prop = "gene_id"
value = value[7:]
if prop == "translation":
# se translation, também guardar o tamanho
dictionary[tag]["length"] = len(value)
if prop in dictionary[tag]:
# se já encontramos esta propriedade
# para este locus_tag
# verificar que é a mesma
current_value = dictionary[tag][prop]
assert value == current_value
else:
# se não, adicionar
dictionary[tag][prop] = value
return dictionary
Neste momento já temos a informação presente no NCBI relativa à nossa zona do genoma.
Falta verificar se esta informação está de acordo com a da tabela presente em https://www.ncbi.nlm.nih.gov/genome/proteins/416?genome_assembly_id=166758.
Reparámos que esta página faz um pedido HTTP a uma API do NCBI, e usamos essa mesma API para retirar as informações que nos interessavam:
def fetch_table():
"""
A página https://www.ncbi.nlm.nih.gov/genome/proteins/416?genome_assembly_id=166758
faz um pedido http ao url_prefix para preencher a tabela.
Este método implementa a paginação que seria feita manualmente
no website, e guarda a informação que queremos num dicionário.
- A coluna 5 corresponde à propriedade db_xref do genbank (gene_id para nós)
- A coluna 6 corresponde à propriedade gene do genbank
- A coluna 7 corresponde à propriedade locus_tag do genbank
- A coluna 8 corresponde à propriedade protein_id do genbank
- A coluna 11 corresponde à propriedade product do genbank
"""
url_prefix = "https://www.ncbi.nlm.nih.gov/genomes/Genome2BE/genome2srv.cgi?action=GetFeatures4Grid&type=Proteins&genome_id=416&genome_assembly_id=166758&gi=&mode=2&is_locus=1&is_locus_tag=1&optcols=1,1,1,0,0&replicons=52840256,NC_002942.5,chr"
end = False
page = 1
page_size = 100
locus_tag_index = 7
column_mapping = {
5: "gene_id",
6: "gene",
8: "protein_id",
11: "product"
}
nested_values = [5, 8]
dictionary = {}
while not end:
# o crawling acaba quando a página que pedi for vazia
try:
url = url_prefix + "&page=" + str(page) + "&pageSize=" + str(page_size)
file_path, _ = urlretrieve(url)
tree = parse_xml(file_path, add_root=True)
rows = tree.findall("TR")
for row in rows:
cols = row.findall("TD")
locus_tag = cols[locus_tag_index].text
# criar um dicionário vazio para esta locus tag
dictionary[locus_tag] = {}
# para cada uma das colunas que nos interessa,
# guardar o valor no dicionário se for
# diferente de -
for index in column_mapping:
prop = column_mapping[index]
if index in nested_values:
value = cols[index].find("a").text
else:
value = cols[index].text
if value != "-":
dictionary[locus_tag][prop] = value
page += 1
end = len(rows) == 0
except Exception as ex:
print(ex)
end = True
return dictionary
Esta função está definida no módulo util.www e utiliza uma outra função parse_xml também presente no mesmo módulo que pode modificar o ficheiro onde se encontra o XML antes de fazer parsing dele, dependendo dos argumentos:
def parse_xml(file_path, add_root=False, start=0, end=0):
"""
Faz parsing de um ficheiro xml e retorna um 'ElementTree'.
"""
if start > 0 or end > 0:
rw.truncate_file(start, end, file_path)
if add_root:
rw.wrap_file("<root>", "</root>", file_path)
fd = open(file_path, "r")
tree = ElementTree.parse(fd)
fd.close()
return tree
Criámos ainda uma função que verifica se os valores presentes na tabela correspondem aos valores provenientes do NCBI. Esta função retorna os valores diferentes para uma posterior análise manual.
def verify(dictionary, table):
"""
Verifica se os valores presentes na tabela correspondem
aos valores presentes no dicionário.
Se encontrar alguma diferença retorna um dicionário
para análise manual dos resultados.
"""
diff = {}
for tag in dictionary:
# para cada uma das locus_tag da nossa zona
if tag in table:
# se estiver na tabela
values = table[tag]
for prop in values:
dict_value = dictionary[tag][prop]
table_value = values[prop]
if not table_value == dict_value:
diff[tag] = {}
diff[tag]["prop"] = prop
diff[tag]["dictionary"] = dict_value
diff[tag]["table"] = table_value
return diff
Durante a análise manual das diferenças concluímos que a informação era a mesma, e, portanto, não houve necessidade de atualizar a que já tínhamos.
Por fim, para sabermos mais sobre a nossa zona do genoma, mapeamos os gene_id retirados do NCBI para identificadores UniProt, adicionando uma nova propriedade uniprot_id.
def retrieve_uniprot_ids(dictionary):
"""
Adicionar ao dictionário uma nova propriedade:
- uniprot_id
"""
gene_id_to_tag = {}
for tag in dictionary:
if "gene_id" in dictionary[tag]:
gene_id = dictionary[tag]["gene_id"]
gene_id_to_tag[gene_id] = tag
gene_ids = list(gene_id_to_tag.keys())
gene_id_to_uniprot_id = www.gene_ids_to_uniprot_ids(gene_ids)
for gene_id in gene_id_to_uniprot_id:
uniprot_id = gene_id_to_uniprot_id[gene_id]
tag = gene_id_to_tag[gene_id]
# adicionar uma nova propriedade ao dicionário: "uniprot_id"
dictionary[tag]["uniprot_id"] = uniprot_id
return dictionary
Todos estes passos podem ser encontrados no ficheiro first.py e resumem-se a:
import util.www as www
start = 270001
end = 505535
record = www.fetch_genbank(start, end)
features = extract_features(record)
dictionary = features_to_dictionary(start, features)
table = www.fetch_table()
diff = verify(dictionary, table)
# análise manual de diff
dictionary = retrieve_uniprot_ids(dictionary)
Na verdade, o que podemos encontrar no ficheiro first.py é:
import util.www as www
import util.rw as rw
start = 270001
end = 505535
ncbi_gb_path = ".ncbi.gb"
ncbi_json_path = ".ncbi.json"
#record = www.fetch_genbank(start, end)
#rw.write_genbank(record, ncbi_gb_path)
record = rw.read_genbank(ncbi_gb_path)
features = extract_features(record)
dictionary = features_to_dictionary(start, features)
rw.write_json(dictionary, ncbi_json_path)
#...