Visualização de dados em Python

Técnicas de visualização de dados podem expressar grandes quantidades de informações de maneira eficiente e auxiliar na detecção de padrões que, sem elas, seriam praticamente imperceptíveis. Essas técnicas utilizam elementos visuais para evidenciar padrões em dados, que conseguimos assimilar melhor devido ao fato de nosso cérebro tender a criar associações harmônicas entre formas e cores. Como dizia aquele velho ditado: “uma imagem vale mais do que mil palavras”.

Em Bioinformática, a visualização de dados surge como um elemento transformador, capaz de representar a exagerada massa de dados produzida pelas recentes evoluções em equipamentos biotecnológicos. Visualização de dados também fornece meios para visão dos processos naturais em pequenas estruturas que compõe a vida, como as interações entre diferentes tipos de moléculas.

Python fornece diversas ferramentas para criação de figuras e gráficos, como por exemplo, histogramas, gráficos de linhas, pontos e frequência. Além disso, Biopython fornece classes e métodos para visualização de dados biológicos, visualização sequências e, até mesmo, de genomas completos.

Nos capítulos anteriores, realizamos o download de um arquivo no formato GenBank com anotações gênicas de um plasmídeo (NC_009934.gbk). Para os exemplos a seguir precisaremos de um arquivo um pouco maior. Para esses exemplos, usaremos o arquivo GenBank com o genoma completo de Acetobacter pasteurianus IFO 3283-12 (NC_017108.gbk). Você pode obtê-lo diretamente no banco de dados GenBank do NCBI. Entretanto disponibilizamos uma cópia do arquivo “NC_017108.gbk” no diretório exemplos do repositório deste livro no GitHub: <https://github.com/dcbmariano/biopython>.

Histogramas

Histogramas são gráficos de distribuição de frequência. Com o pacote Bio.Seq podemos criar uma lista com tamanhos de sequências codificantes presentes no nosso arquivo exemplo (NC_017108.gbk) e, em seguida, usar a módulo pylab para plotá-las.

from Bio import SeqIO
import pylab
 
# Recebendo arquivo GBK
exemplo = SeqIO.read("NC_017108.gbk", "genbank")
tamanhos = []
 
# Obtendo quantidade de aminoacidos em sequencias
for i in exemplo.features:
	if i.type == 'CDS':
		tamanhos.append(len( \
i.qualifiers['translation'][0]))
 
# Gerando a figura
pylab.hist(tamanhos, bins=20)
pylab.title("Histograma - frequencia de sequencias")
pylab.xlabel("Tamanho da sequencia (bp)")
pylab.ylabel("Quantidade")
pylab.show()

Inicialmente foi realizada a importação dos módulos necessários e a leitura do arquivo GenBank. Observe que uma lista chamada tamanhos foi declarada em branco ([]). Essa lista foi utilizada para armazenar os tamanhos das strings presentes em i.qualifiers[‘translation’][0]. A partir desse ponto, alguns módulos do pylab foram utilizadas para gerar a figura:

• hist( ): recebe a lista todos os tamanhos. O parâmetro bins define a quantidade de barras do histograma;
• title( ): define o título do gráfico;
• xlabel( ): define o título do eixo x;
• ylabel( ): define o título do eixo y;
• show( ): gera a figura. O programa continuará executando enquanto a figura não for fechada.

O histograma nos mostra que em nosso arquivo exemplo há uma quantidade maior de genes que codificam proteínas com tamanho de aproximadamente 200 aminoácidos. Ainda nos mostra a existência de genes que codificam proteínas com tamanho próximo a 2.000 aminoácidos. Como seriam necessários três nucleotídeos para codificar um aminoácido, isso indica a existência de regiões codificadores maiores que 6.000 pares de base.

Gráfico de pizza

No exemplo a seguir utilizaremos o método pie( ) para criar um gráfico de pizza (ou gráfico de torta) sobre o conteúdo GC. O conteúdo GC é a quantidade de bases nitrogenadas no DNA que são ou guanina ou citosina.

from Bio import SeqIO
from Bio.SeqUtils import GC
import pylab
 
for i in SeqIO.parse("NC_017108.gbk", "genbank"):
	gc = GC(i.seq)
 
at = 100-gc
 
pylab.pie([gc,at])
pylab.title("Conteudo GC:")
pylab.xlabel("GC: %0.1f porcento\nAT: %0.1f porcento" \
% (gc,at))
pylab.show()

Nesse exemplo, inicialmente foi importado o módulo GC do pacote Bio.SeqUtils que permitirá a análise de conteúdo GC. O conteúdo GC é extraído através do método GC( ). O gráfico de pizza é criado através do método pie( ), que recebe como entrada uma lista com as variáveis gc (conteúdo GC) e at (criada pela subtração do valor máximo possível, 100 %, pelo conteúdo GC).

Entretanto, a visualização gerada não nos apresenta um resultado relevante, uma vez que a simples informação obtida pelo valor do conteúdo GC (53%) já seria suficiente compreender o necessário sobre o conteúdo GC total. Podemos obter visualizações mais elegantes e informativas do conteúdo GC utilizando gráficos de linhas.

Gráfico de linhas

Utilizando a sequência de genoma completo do arquivo, vamos dividi-lo em fragmentos de mil pares de base e plotar o conteúdo GC de cada fragmento em um gráfico de linhas.

from Bio import SeqIO
from Bio.SeqUtils import GC
import pylab
 
# Parte 1 - Abrindo arquivo GBK
for i in SeqIO.parse("NC_017108.gbk", "genbank"):
	seq = str(i.seq)
 
# Parte 2 - Variaveis importantes
tamanho = len(seq)
fragmentos = int(tamanho / 1000)
gc = []
 
# Parte 3 - Armazenando conteudo GC
for i in range(fragmentos):
	j = i * 1000
	k = j + 999
	gc_atual = GC(seq[j:k])
	gc.append(gc_atual)
	print(i,": ",j,"-",k,"- GC =",gc_atual,"%" )
 
# Parte 4 - Adicionando o ultimo elemento
resto = tamanho % 1000
j = (i+1) * 1000
k = j + resto
gc_ultimo = GC(seq[j:k])
gc.append(gc_ultimo)
print(i+1,": ",j,"-",k,"- GC =",gc_ultimo,"%" )
 
# Parte 5 - Imprimindo grafico
pylab.plot(gc)
pylab.title("Conteudo GC\n%i fragmentos de 1000 pb \
variando de %0.1f%% a %0.1f%%" % (len(gc),min(gc),max(gc)))
pylab.xlabel("Tamanho (Kb)")
pylab.ylabel("GC%")
pylab.show()

Para uma melhor explicação, dividimos o código em cinco partes:

Parte 1: após a importação dos módulos necessários, o arquivo GenBank é analisado, e o genoma completo é salvo como string na variável seq.

Parte 2: o tamanho do genoma e a quantidade de fragmentos de exatos mil pares de base são determinados. Além disso, uma lista vazia chamada gc é criada. Será nessa lista que os valores de conteúdo GC serão armazenados.

Parte 3: a partir desse ponto, foi realizada a determinação do conteúdo GC por fragmento. Observe que para isso foram criadas três variáveis: i, j e k. A variável i percorre todos os fragmentos. A variável j recebe a posição inicial do fragmento, enquanto k recebe a final. Os valores são armazenados e impressos na tela:

0 :  0 - 999 - GC = 52.5525525526 %
1 :  1000 - 1999 - GC = 49.049049049 %
2 :  2000 - 2999 - GC = 48.4484484484 %
3 :  3000 - 3999 - GC = 49.9499499499 %
...

Parte 4: como determinamos que o laço só irá percorrer fragmentos de tamanho 1000, o último elemento, por provavelmente ter tamanho inferior a 1000, não será percorrido. Por isso usamos o operador módulo (%) para determinar seu tamanho e assim armazená-lo na lista.

2904 :  2904000 - 2904624 - GC = 53.2051282051 %

Parte 5: por fim, o método pylab.plot(gc) gera o gráfico, os métodos title, ylabel e xlabel formatam a figura, e show exibe o gráfico na tela.

A grande variação do conteúdo GC nos dados dificulta a visualização do resultado. Entretanto, você pode alterar o tamanho de cada fragmento e obter uma visualização com variação de linhas mais suave.

Na figura acima, cada fragmento equivale a 10.000 pb.

Semelhanças entre sequências através de gráficos de pontos

Neste exemplo, iremos avaliar semelhanças entre duas sequências utilizando gráficos de pontos. Crie um arquivo FASTA (exemplo_pontos.fasta):

>beta-glucosidase_parcial_mRNA
CGTGTGCTTGCGCTACTGCTTGGTCTCACTGTGTGAAC
AACTGCTGCCGCCGCCTGCCGGCGGTGTGCTC
>beta-glucosidase_parcial_mRNA2
CGTGTGCTTGCGCTACTGCTTGGTCTCACTGTGTGACC
AACTGCTGCCGCCGCCTGCCGGCGGTGTGCTC

Neste exemplo, recomendado pela documentação oficial do Biopython, iremos ler as duas sequências FASTA e percorrer as sequências usando janela de tamanho 11.

from Bio import SeqIO
import pylab
 
handle = open("exemplo_pontos.fasta")
record_iterator = SeqIO.parse(handle, "fasta")
rec_one = next(record_iterator)
rec_two = next(record_iterator)
handle.close()
 
# Definindo tamanho de janela
window = 11
dict_one = {}
dict_two = {}
 
# Comparacoes

for (seq, section_dict) in \
[(str(rec_one.seq).upper(), dict_one), \
(str(rec_two.seq).upper(), dict_two)]:
    for i in range(len(seq)-window):
        section = seq[i:i+window]
        try:
            section_dict[section].append(i)
        except KeyError:
            section_dict[section] = [i]
 
matches = set(dict_one).intersection(dict_two)
print("%i k-mers identicos" % len(matches)) )
 
# Determinando posicoes x e y
x = []
y = []
for section in matches:
    for i in dict_one[section]:
        for j in dict_two[section]:
            x.append(i)
            y.append(j)
 
pylab.cla() 
pylab.gray()
pylab.scatter(x,y)
pylab.xlim(0, len(rec_one)-window)
pylab.ylim(0, len(rec_two)-window)
pylab.xlabel("%s (tamanho %i pb)" \
% (rec_one.id, len(rec_one)))
pylab.ylabel("%s (tamanho %i pb)" \
% (rec_two.id, len(rec_two)))
pylab.title("Grafico de pontos com janela \
de %i pb\n(sem mismatches)" % window)
pylab.show()

Nesse exemplo realizamos a análise do arquivo Multi-FASTA. Lemos duas sequências (rec_one e rec_two) e realizamos a comparação entre elas usando uma janela (window) de tamanho 11, ou seja, comparamos de 11 em 11 caracteres.

O trecho em que são realizadas as comparações entre sequências é bastante complexo. Destacamos que a variável section recebe os fragmentos da sequência de 11 caracteres e os dicionários dict_one e dict_two que recebem a posição da janela e as sequências. Observe o uso das funções try, para tentar executar um comando, e except, caso algum comando presente em try falhe.

Em seguida, o script determina as posições x e y em que serão impressos os pontos e os métodos da biblioteca pylab são utilizados para gerar a figura.

Uma das falhas graves dessa análise é a não detecção de semelhanças em sequências em reverso complementar. Os autores da documentação do Biopython relatam a necessidade de aperfeiçoar esse script inserindo esse tipo de análise, e incentivam sua moficação e divulgação por parte de usuários a fim de contribuir com a comunidade de desenvolvedores do Biopython.

Entretanto, ressaltamos que existem métodos mais eficientes para comparação entre sequências e até mesmo genomas completos, como por exemplo, utilizando o pacote Bio.Graphics. Entretanto o pacote Bio.Graphics requer a instalação de algumas bibliotecas de terceiros. Abordaremos isso no próximo tópico.

Instalando bibliotecas de terceiros para visualização de genomas

Biopython ainda fornece o pacote Bio.Graphics, que depende da biblioteca ReportLab. Para instalá-la, acesse o Website da ReportLab <http://www.reportlab.org/> e registre uma conta gratuitamente. Depois faça o download da última versão da biblioteca em: <https://www.reportlab.com/pypi/simple/reportlab/>.

No Windows realize a instalação do arquivo executável. No Linux e MacOS utilize os comandos abaixo para descompactar e instalar:

tar -zxvf reportlab-INSIRA_A_VERSAO_AQUI.tar.gz
sudo python setup.py install

Instale também a biblioteca Python Imaging Library (PIL). Para instalá-la acesse: <http://www.pythonware.com/products/pil/> e faça o download da última versão compatível com seu sistema operacional.

Visualização de genomas

Neste primeiro exemplo, usaremos como exemplo o arquivo “NC_009934.gbk” (confira no diretório exemplos). Vamos imprimir na tela o genoma de maneira circular. Veja:

from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics import GenomeDiagram
from Bio import SeqIO
 
# Recebendo o arquivo
arquivo = SeqIO.read("NC_009934.gbk", "genbank")
 
# Configuracoes basicas
gd_diagram = GenomeDiagram.Diagram("Plasmideo")
gd_track_for_features = gd_diagram.new_track(1, \
name="Anotacoes")
gd_feature_set = gd_track_for_features.new_set()
 
# Imprime as features de maneira intercalada
for feature in arquivo.features:
    if feature.type != "gene":
        # Despreza a feature que nao sao genes
        continue
    if len(gd_feature_set) % 2 == 0:
        color = colors.HexColor('#79B134')
    else:
        color = colors.HexColor('#8DE91D')
    gd_feature_set.add_feature(feature, color=color, \
label=True)
 
# Desenhando

gd_diagram.draw(format="circular", circular=True, \
pagesize=(20*cm,20*cm),start=0, end=len(arquivo), \
circle_core=0.7)
 
# Salvand em 2 formatos: PDF e PNG
gd_diagram.write("plasmideo.pdf", "PDF")
gd_diagram.write("plasmideo.png", "PNG")

Ao executar esse código dois arquivos serão produzidos: um no formato PNG e outro no formato PDF.

Agora entenda um pouco do que foi feito:

1. importamos os pacotes e módulos necessários, como módulo colors do pacote reportlab.lib, responsável pelas funções relacionadas a cores, cm do reportlab.lib.units, responsável por identificar os valores de medidas de páginas, o GenomeDiagram do Bio.Graphics, responsável pelas funções de impressão de genomas, e o SeqIO;
2. recebemos o arquivo no formato GenBank;
3. no trecho de configurações básicas, declaramos um objeto da classe GenomeDiagram, que criará um novo diagrama chamado de “plasmídeo”, adicionamos uma nova faixa e um conjunto de recursos;
4. em seguida, usamos um laço para percorrer todas as features. Criamos um trecho em que imprimimos somente os genes com cores intercaladas. Você pode brincar com as cores trocando o valor do código hexadecimal, por exemplo, usamos “color = colors.HexColor(‘#8DE91D’)”, que é uma variação do tom verde, você poderia alterar o valor #8DE91D para #990000 que é uma variação de vermelho, ou poderia simplesmente escrever o nome da cor, como por exemplo “color = colors.red” (você pode aperfeiçoar as cores pesquisando mais sobre o código de cores hexadecimal RGB);
5. o método draw( ) realiza o desenho recebendo como parâmetros o formato (format=”circular”, circular=True), o tamanho da página (pagesize=(20*cm,20*cm)), a posição de início e fim (start=0, end=len(arquivo)) e a distância do centro do círculos (circle_core=0.7); e
6. por fim, salvamos em dois formatos: PNG (formato de imagens simples) e PDF (formato vetorizado).

Você ainda pode alterar o modo de exibição dos genes e usar setas para identificar o sentido de transcrição. Para realizar isso é necessário inserir o parâmetro sigil=”ARROW” dentro do método gd_feature_set.add_feature( ).

gd_feature_set.add_feature(feature, color=color, \
label=True,sigil="ARROW")

Podemos ainda, modificar nosso código para que ele exiba o resultado de maneira linear. Para isso é necessário apenas alterar o trecho onde o desenho é construído. Modifique o trecho:

# Desenhando 

gd_diagram.draw(format="circular", circular=True, \
pagesize=(20*cm,20*cm),start=0, end=len(arquivo),  \
circle_core=0.7)

Por este trecho:

# Desenhando - figura linear

gd_diagram.draw(format="linear", orientation="landscape", \
pagesize='A4',fragments=1, start=0, end=len(arquivo))

Nesse exemplo, as principais alterações foram nos parâmetros: format (“circular” para “linear”), orientation (foi inserido o valor landscape para que o papel ficasse na horizontal), pagesize (foi inserido o valor “A4”, que corresponde a uma folha de tamanho 21 cm por 29,7 cm) e fragment (que corresponde ao número de linhas a qual o genoma será exibido). Executando o código, obteremos uma figura assim:

Visualizando marcações em genomas

Nos exemplos anteriores, visualizamos genes presentes em genomas pequenos, como no caso, de plasmídeos. Se quiséssemos visualizar grandes genomas, os textos que identificam os genes se iriam sobrepor, impossibilitando a leitura e identificação de genes. Veja por exemplo, a figura circular gerada pelo script descrito anteriormente e com o arquivo “NC_017108.gbk”.

Podemos compreender que é importante identificar apenas elementos relevantes para determinadas análises, como por exemplo, identificar RNAs estruturais.

No exemplo a seguir buscaremos identificar RNAs estruturais e demonstrá-los na figura a partir de dois tipos de marcadores: (i) a cor azul para genes que codificam RNAs ribossomais (rRNA); e (ii) a cor verde genes que codificam RNAs transportadores (tRNA).

Primeiro, vamos visualizar apenas os tRNAs. Veja como ficaria nosso código:

from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics import GenomeDiagram
from Bio import SeqIO
 
# Recebendo o arquivo
arquivo = SeqIO.read("NC_017108.gbk", "genbank")
 
# Configuracoes basicas
gd_diagram = GenomeDiagram.Diagram("Marcadores")
gd_track_for_features = gd_diagram.new_track(1, \
name="Anotacoes")
gd_feature_set = gd_track_for_features.new_set()
 
for feature in arquivo.features:
 
	if feature.type == "tRNA":
		color = colors.green
		feature.qualifiers['locus_tag'] = \
feature.qualifiers['product']
		gd_feature_set.add_feature(feature, \
color=color, label=True)
 
	else:
		continue
 
# Desenhando
gd_diagram.draw(format="circular", circular=True, \
pagesize=(20*cm,20*cm),start=0, end=len(arquivo), \
circle_core=0.7)
 
# Salvand em 2 formatos: PDF e PNG
gd_diagram.write("trna.pdf", "PDF")
gd_diagram.write("trna.png", "PNG")

Observe que validamos se a feature pertence a um tRNA em “if feature.type == “tRNA”:”. Após essa validação, declaramos a cor como verde e inserimos um comando que você deve ter considerado um pouco curioso: “feature.qualifiers[‘locus_tag’] = feature.qualifiers[‘product’]”. Declaramos que o valor presente em product deve ser transferido ao campo locus_tag. Isso é necessário neste caso, pois na hora de gerar a figura, a biblioteca utiliza o valor presente no campo locus_tag como legenda. Por uma questão de estética achamos melhor imprimir o valor presente no campo product.

Observe que a figura gerada ao final apenas apresenta marcações nas posições onde estão os tRNAs.

Para visualizar a localização de rRNAs, altere os trechos:

if feature.type == "tRNA":
	color = colors.green
	feature.qualifiers['locus_tag'] = \
feature.qualifiers['product']
	gd_feature_set.add_feature(feature, \
color=color, label=True)

Por:

if feature.type == "rRNA":
	color = colors.blue
	gd_feature_set.add_feature(feature, color=color, \
label=True)

Altere também o nome dos arquivos salvos (mude para rrna.pdf e rrna.png).

Observe que agora somente são exibidas as posições de rRNAS. Podemos ainda, unir os resultados de tRNAs e rRNAs, e inserir marcações coloridas para detectar a posição dos genes. Veja:

from reportlab.lib import colors
from Bio.Graphics import GenomeDiagram
from Bio import SeqIO
 
# Recebendo o arquivo
arquivo = SeqIO.read("NC_017108.gbk", "genbank")
 
# Configuracoes basicas
gd_diagram = GenomeDiagram.Diagram("Marcadores")
gd_track_for_features = gd_diagram.new_track(1, \
name="Anotacoes")
gd_feature_set = gd_track_for_features.new_set()
 
# Imprime as features de maneira intercalada
for feature in arquivo.features:
 
	if feature.type == "rRNA":
		color = colors.blue
		gd_feature_set.add_feature(feature, \
color=color, label=True)
 
	elif feature.type == "tRNA":
		color = colors.green
		feature.qualifiers['locus_tag'] = \
feature.qualifiers['product']
		gd_feature_set.add_feature(feature, \
color=color, label=True)
 
	elif feature.type == "gene":
 
		if len(gd_feature_set) % 2 == 0:
			color = colors.HexColor('#BBBBBB')
		else:
			color = colors.HexColor('#DDDDDD')
		
		gd_feature_set.add_feature(feature, \
color=color, label=False)
 
	else:
		continue
 
# Desenhando
gd_diagram.draw(format="circular", circular=True,\
 pagesize=(20*cm,20*cm),start=0, end=len(arquivo), \
circle_core=0.7)
 
# Salvand em 2 formatos: PDF e PNG
gd_diagram.write("completo.pdf", "PDF")
gd_diagram.write("completo.png", "PNG")

Observe que foi realizado um comando condicional para detectar cada tipo de estrutura. Também intercalamos as cores de impressão de genes. Para isso utilizamos dois tons diferentes de cinza: (i) “#BBBBBB” e (ii) “#DDDDDD”.

Sintenia entre genomas

Sintenia pode ser definida como: quando blocos gênicos são conservados em cromossomos de dois ou mais organismos diferentes quando comparados entre si. Para realizar análises de sintenia entre genomas pode-se utilizar a ferramenta BLAST para comparar as sequências entre eles.

No exemplo a seguir, utilizaremos os conhecimentos aprendidos até aqui sobre manipulação de arquivos e comparações de sequências utilizando BLAST através do Biopython, e juntá-los aos conhecimentos ensinados nesse capítulo sobre visualização de dados.

Neste exemplo, receberemos dois arquivos de genomas completos, cujos nomes devem ser inseridos pelo usuário, detectaremos as similaridades entre eles e geraremos um gráfico de sintenia.

Vamos exibir o código fonte e realizar estudo de caso. A seguir, iremos explicar os pontos mais importantes do código.

from Bio import SeqIO
from Bio.Graphics import GenomeDiagram
from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics.GenomeDiagram import CrossLink
from reportlab.lib import colors
from Bio.Blast.Applications import *
 
import sys
 
# Recebendo os dados pela chamada do programa
A = sys.argv[1]
B = sys.argv[2]
 
A1 = SeqIO.read(A,"genbank")
B1 = SeqIO.read(B,"genbank")
 
# BLAST
SeqIO.convert(A, "genbank", A+".fasta", "fasta")
SeqIO.convert(B, "genbank", B+".fasta", "fasta")


comando_blastn = NcbiblastnCommandline( \
query=A+".fasta", subject=B+".fasta", \
outfmt="'6 qstart qend sstart send pident'",\
out="blast_"+A+"_"+B+".txt")
stdout, stderr = comando_blastn()
blast = open("blast_"+A+"_"+B+".txt")
 
# Iniciando a figura
name = A+"_"+B
gd = GenomeDiagram.Diagram(name)
 
gA = gd.new_track(1,name="A",height=0.5, \
start=0,end=len(A1))
gA1 = gA.new_set()
gB = gd.new_track(3,name="B",height=0.5, \
start=0,end=len(B1))
gB1 = gB.new_set()
 
# Cores CDSs - intercalado
c1 = "#79B134"
c2 = "#8DE91D"
 
# Colore um quadrado para cada CDS do arquivo A
cont = 1
for i in A1.features:
 
	if i.type == "CDS":
 
		if cont % 2 == 1:
			color_atual = c1
		else:
			color_atual = c2
		gA1.add_feature(i, label=False, \
label_position="start",color=color_atual)
 
		cont += 1
 
	if i.type == "rRNA":
		color_atual = colors.blue
		gA1.add_feature(i, label=False,\
label_position="start",color=color_atual)
 
# Colore um quadrado para cada CDS do arquivo B
cont = 1
for i in B1.features:
 
	if i.type == "CDS":
		if cont % 2 == 1:
			color_atual = c1
		else:
			color_atual = c2
		gB1.add_feature(i, label=False,\
label_position="start",color=color_atual)
		cont += 1
 
	if i.type == "rRNA":
		color_atual = colors.blue
		gB1.add_feature(i, label=False,\
label_position="start",color=color_atual)
 
# Marca na figura os trechos sintenicos
for b in blast:
	qstart = int(b.split("\t")[0])
	qend = int(b.split("\t")[1])
	sstart = int(b.split("\t")[2])
	send = int(b.split("\t")[3])
	identidade = (float(b.split("\t")[4])*0.8)/100
 
	# Detectando inversoes
	qinv = qend - qstart
	sinv = send - sstart
 
	if (qinv > 0 and sinv > 0) or \
(qinv < 0 and sinv < 0):
		cor = colors.Color\
(1,.341176,.341176,identidade)
	else:
		cor = colors.firebrick
		
	gd.cross_track_links.append(CrossLink((gA, \
qstart, qend),(gB, sstart, send),color=cor))
 
gd.draw(format="linear", pagesize=(8*cm,29.7*cm), \
fragments=1)
 
gd.write(name + ".pdf", "PDF")

Para testar o script desenvolvido, vamos realizar um estudo de caso com arquivos GenBank de dois genomas completos: Escherichia coli 042 (NC_008253.gbk) e Escherichia coli 536 (NC_017626.gbk).

Faça o download dos arquivos (pelo GenBank ou pelo diretório exemplos) e do script (diretório scripts do GitHub, arquivo c13_s13.py). Execute o script com a linha de comando:

python c13_s13.py NC_017626.gbk NC_008253.gbk

Como resultado você obterá dois arquivos:

1. NC_017626.gbk_NC_008253.gbk.txt
2. NC_017626.gbk_NC_008253.gbk.pdf

O primeiro é um um arquivo de texto com os resultados de BLAST. O segundo é um arquivo PDF com um gráfico de sintenia entre os dois genomas.

Na parte de cima do gráfico de sintenia é impresso o primeiro genoma enviado como argumento. Os genes são impressos em tons intercalados de verde e, na cor azul, os RNAs ribossomais. Na parte de baixo, também em tons intercalados de verde, está impresso o segundo genoma. As regiões sintenicas entre os dois genomas são interligadas por segmentos de reta de cor vermelha. Segmento na cor vermelha claro indicam regiões sintênicas comuns, enquanto a cor vermelha escura indica inversões gênicas. No exemplo demonstrado, podemos perceber que ambos os genomas são altamente sintênicos.

Vamos agora analisar os principais pontos do código fonte. Vamos começar na parte em que são importados módulos e pacotes:

from Bio.Graphics.GenomeDiagram import CrossLink

O módulo Crosslink será responsável por imprimir os segmentos de reta que interligam um genoma a outro. Na parte de criação do diagrama destacamos o trecho em que são criadas as faixas lineares em que serão impressos os genomas.

gA = gd.new_track(1,name="A",height=0.5, \
start=0,end=len(A1))

gB = gd.new_track(3,name="B",height=0.5, \
start=0,end=len(B1))

As variáveis gA e gB armazenam os objetos que registram as faixas a qual serão armazenados os dados de ambos os genomas. Observe que o método new_track recebe atributos como height, que define a altura da faixa impressa que representam os genomas, start e end, que definem início e fim das faixas (necessário devido ao fato dos genomas possuírem tamanhos diferentes).

# Marca na figura os trechos sintenicos
for b in blast:
	qstart = int(b.split("\t")[0])
	qend = int(b.split("\t")[1])
	sstart = int(b.split("\t")[2])
	send = int(b.split("\t")[3])
	identidade = (float(b.split("\t")[4])*0.8)/100

Esse trecho percorre a lista que armazena os resultados de BLAST e extrai cinco informações principais: (i) qstart; (ii) qend; (iii) sstart; (iv) send; e a (v) identidade. Crosslink fará marcações da posição qstart a qend e as interligará a outras marcações da posição sstart a send.

Ainda utilizamos operações matemática para detectar inversões que, se existirem, as marcações serão impressas na cor firebrick (um vermelho mais forte). As interligações são inseridas no gráfico no trecho:

gd.cross_track_links.append(CrossLink((gA, \
qstart, qend),(gB, sstart, send),color=cor))

A identidade é utilizada apenas para definir o nível de transparência da cor impressa (efeito alfa). Observe que multiplicamos por um fator de 0.8 para forçar uma transparência mínima de 20% e dividimos por 100, pois o resultado de BLAST é dado em porcentagem.

Podemos ainda modificar nosso código para que apenas imprima trechos em que a identidade é igual a 100 %. Para isso vamos apagar o “fator de transparência” e inserir um comando condicional antes da impressão dos segmentos de retas que interligam os genomas pelo Crosslink.

Altere a linha:

identidade = (float(b.split("\t")[4])*0.8)/100

Para:

identidade = (float(b.split("\t")[4]))/100

E altere o trecho antes da linha:

gd.cross_track_links.append(CrossLink((gA, \
qstart, qend),(gB, sstart, send),color=cor))

Inserindo o comando condicional:

if identidade >= 1:
	gd.cross_track_links.append(CrossLink((gA, \
qstart, qend),(gB, sstart, send),color=cor))

Ao executar o novo código, você obterá como resultado uma figura com uma quantidade maior de trechos em branco, entretanto será possível identificar uma maior quantidade de informações, como por exemplo, que o começo do genoma exibido acima possui uma maior similaridade com o final do genoma abaixo, e vice-versa. Além de que no centro há um trecho com uma sintenia menor.

Você pode alterar o código e inserir outros recursos, como nomes de genes, alterar as cores utilizadas ou até mesmo inserir outras marcações para RNAs estruturais. Também é possível obter visualizações mais próximas de genes, alterando os campos start e end do método new_track.

Biopython fornece ainda muitos outros métodos para visualização de dados. Você pode saber mais acessando a documentação oficial. Você pode ainda utilizar a função help( ) para saber mais sobre os métodos apresentados neste capítulo.

Quer aprender mais? Conheça nossos cursos profissionalizantes à partir de R$19,99: