Outras coisas interessantes que se pode fazer com Biopython

Neste capítulo serão apresentadas alguns exemplos de tarefas interessantes que podem ser implementadas com Biopython. Alguns exemplos estão presentes na documentação oficial (e podem ser visualizados no capítulo “Cookbook – Cool things to do with it”), entretanto outros foram desenvolvidos com base no que foi apresentado nos capítulos anteriores. Ressaltamos ainda que neste capítulo não iremos focar em explicações detalhadas sobre cada ponto do código, mas sim em demonstrações de códigos interessantes que podem ser desenvolvidos com Biopython.

Ordenando sequências por tamanho

Vamos criar um script que recebe um arquivo Multi-FASTA e ordena as sequências pelo tamanho.

Primeiramente, crie um arquivo de sequências Multi-FASTA ou faça download de um, como por exemplo, as sequências de nucleotídeos codificantes de Acetobacter pasteurianus (NC_017108.ffn). Agora vamos utilizar o método SeqIO.parse( ) para ler o arquivo, o método sort( ) para ordenar e o método SeqIO.write( ) para gravar as sequências ordenadas por tamanho em um arquivo FASTA.

from Bio import SeqIO
 
# Lendo as sequencias
sequencias = list(SeqIO.parse("NC_017108.ffn","fasta"))
 
# Realizando o ordenamento
sequencias.sort(cmp=lambda x,y: cmp(len(x),len(y)))
 
# Gravando o resultado 
if SeqIO.write(sequencias, "sequencias_ordenadas.fasta", \
"fasta"):
	print("Sequencias ordenadas com sucesso." )
else:
	print("Um erro ocorreu." )

Calculando N50 usando Biopython

O N50 é uma métrica utilizada para avaliar montagens de genomas. N50 constitui no tamanho da sequência a qual está presente o par de base que corresponde a 50% do total de bases em um arquivo Multi-FASTA ordenado do menor para o maior. No script abaixo usaremos os comandos aprendidos anteriormente para calcular o valor de N50 de um arquivo Multi-FASTA.

# Calculo de N50
 
# Declaracoes iniciais
from Bio import SeqIO
import sys
 
sum_tamanhos = 0
sumcontig = 0
todos_tamanhos = list()
 
# Recebemos o arquivo como parametro
arquivo = sys.argv[1]
 
# Laco le todos os tamanhos 
for i in SeqIO.parse(arquivo,"fasta"):
	tam_contig = len(i.seq)
	todos_tamanhos.append(tam_contig)
	sumcontig = sumcontig + tam_contig
 
# Ordenamos todos os tamanhos
todos_tamanhos.sort()
 
# Calcula metade da soma das sequencias
v50 = sumcontig/2
 
# Laco para detectar N50
for tam in todos_tamanhos:
	sum_tamanhos = sum_tamanhos + tam
	if sum_tamanhos > v50:
		n50 = tam
		break
 
print("N50: ",n50)

No script demonstrado acima, inicialmente foi declarado os módulos que serão utilizados. O nome do arquivo a ser analisado é recebido como parâmetro. Um laço varre o arquivo Multi-FASTA, calcula o tamanho de todas as sequências, os armazena numa lista denominada “todos_tamanhos” e a soma do tamanho de todas as sequências é armazenada em uma variável chamada sumcontig. Em seguida, a lista é ordenada através do método sort( ), então é calculado o valor da metade da soma de todas as sequências. Um laço varre a lista e o script usa um comando condicional pra detectar o tamanho da maior sequência da qual a soma seja inferior a metade da soma das sequências, ou seja, o valor de N50. Por fim, o comando print imprime o resultado na tela.

Contando a quantidade de leituras em um arquivo FASTQ

O formato FASTQ incorpora sequências seguidas por um fator de qualidade, dado pelo algoritmo de PHRED, que indica em uma base logaritma a probabilidade de erro uma base obtida por sequenciamento.

Você pode encontrar exemplos de arquivos FASTQ no banco de dados SRA do NCBI: <http://www.ncbi.nlm.nih.gov/sra>. Como exemplo, faça o download do arquivo “SRR020192.fastq.gz” e extraia o arquivo FASTQ (no Linux e MacOS, apenas abra o arquivo, ou no Windows utilize o software Winrar).

Agora, vamos analisar quantas sequências existem no arquivo FASTQ usando uma variável contadora e a função SeqIO.parse( ).

from Bio import SeqIO
 
# Declarando variavel contadora i
i = 0
 
# Contando registros
for seq in SeqIO.parse("SRR020192.fastq", "fastq"):
    i += 1
 
print("Foram detectadas %i leituras." % i)
 
# Foram detectadas 41892 leituras.

Visualizando a qualidade de bases em arquivos FASTQ

Utilizando os módulos Pylab e SeqIO podemos analisar arquivos do tipo FASTQ e gerar gráficos de qualidade dados pelo algoritmo de PHRED. Neste exemplo, extrairemos a média da qualidade PHRED por posição de par de base em uma leitura, além dos maiores e menores valores de qualidade. Em seguida, iremos gerar um gráfico de qualidade.

from Bio import SeqIO
import pylab
 
arquivo = open("SRR020192.fastq")
parse = list(SeqIO.parse(arquivo,"fastq"))
tam_seq = []
max_seqs = []
min_seqs = []
media = []
 
for i in parse:
	tam_seq.append(len(i.seq))
 
max_seq = max(tam_seq)
min_seq = min(tam_seq)
total_seq = len(parse)
 
for j in range(max_seq):
	aux = 0
	mas = 0
	mis = max_seq
	for i in parse:
		if(len(i.seq) < max_seq):
			i.letter_annotations\
["phred_quality"].append(0)
		aux = aux + \
i.letter_annotations["phred_quality"][j]
 
		if mas < i.letter_annotations\
["phred_quality"][j]:
			mas = i.letter_annotations\
["phred_quality"][j]
 
		if mis > i.letter_annotations\
["phred_quality"][j]:
			mis = i.letter_annotations\
["phred_quality"][j]
 
	media.append(int(aux/total_seq))
	max_seqs.append(mas)
	min_seqs.append(mis)
 
pylab.plot(media,'-k',label="Media", linewidth=2)
pylab.plot(max_seqs,'-g',label="Max")
pylab.plot(min_seqs,'-r',label="Min")
pylab.grid(color='gray', linestyle=':', linewidth=1)
pylab.legend()
pylab.ylim(0,45)
pylab.ylabel("PHRED")
pylab.xlabel("Leitura (pb)")
pylab.savefig("grafico_qualidade.png")
 
pylab.show()

No gráfico de qulidade PHRED, bases com qualidade abaixo de 20 possuem uma alta provabilidade de erro (maior que 1 chance em 100). Leituras com qualidade abaixo de 10 possuem provabilidade de erro de 1 em 10.

Analisando o gráfico gerado pela análise do arquivo “SRR020192.fastq” podemos perceber uma perda de qualidade em média após a posição 125 entre todas as leituras. Essa perda de qualidade pode ainda estar relacionada a presença de uma maior quantidade de leituras menores do que 125 pares de bases. Após a posição 250, a qualidade média de leituras cai a valores próximos a zero.

Criando um filtro de qualidade para arquivos FASTQ

Utilizando o mesmo arquivo do exemplo anterior, vamos agora criar um filtro de qualidade para remover sequências que tenham qualidade abaixo de PHRED 20 (99% de acurácia).

from Bio import SeqIO
 
boas_leituras = (rec for rec in \
              SeqIO.parse("SRR020192.fastq", "fastq") \
              if min(rec.letter_annotations \
["phred_quality"]) >= 20)
 
# Contador i
i = SeqIO.write(boas_leituras, "boa_qualidade.fastq", \
"fastq")
 
print("Foram salvas %i leituras com qualidade \
maior ou igual a 20." % i)
 
# Foram salvas 20050 leituras com qualidade >= a 20.

Removendo primers de arquivos com dados brutos de sequenciamento

No exemplo a seguir faremos a remoção da sequência primer do arquivo FASTQ. Neste exemplo, a sequência “GATGACGGTGT” se repete no início de todas as leituras, logo é necessário removê-la.

Faremos isso em quatro etapas: (i) primeiro detectaremos quais sequências possuem primers; (ii) em seguida, verificaremos se existem sequências sem primers; (iii) criaremos um novo arquivo FASTQ com as sequências cortadas; e (iv) uniremos o resultado.

from Bio import SeqIO
 
# Detectando sequencias que tem primer
primer_reads = (rec for rec in \
                SeqIO.parse("SRR020192.fastq", "fastq") \
                if rec.seq.startswith("GATGACGGTGT"))
SeqIO.write(primer_reads, "com_primer.fastq", "fastq")
 
# Detectando sequencias sem primer
 
trimmed_primer_reads = (rec[11:] for rec in \
                        SeqIO.parse("SRR020192.fastq", \
                        "fastq") \
                        if rec.seq.startswith(\
                        "GATGACGGTGT"))\
SeqIO.write(trimmed_primer_reads, \
"com_primer_cortado.fastq", "fastq")
 
# Criando um novo FASTQ
def trim_primer(record, primer):
    if record.seq.startswith(primer):
        return record[len(primer):]
    else:
        return record
 
trimmed_reads = (trim_primer(record, "GATGACGGTGT") \
                 for record in \
                 SeqIO.parse("SRR020192.fastq", "fastq"))
SeqIO.write(trimmed_reads, "cortado.fastq", "fastq")
 
# Unindo os resultados
def trim_primers(records, primer):
	# Remove primers no inicio das leituras
    len_primer = len(primer) 
 
    for record in records:
        if record.seq.startswith(primer):
            yield record[len_primer:]
        else:
            yield record
 
original_reads = SeqIO.parse("SRR020192.fastq", "fastq")
trimmed_reads = trim_primers(original_reads, "GATGACGGTGT")
i = SeqIO.write(trimmed_reads, "cortado.fastq", "fastq")
 
print("Foram salvas %i leituras" % i)

Convertendo arquivos FASTA e QUAL em um arquivo FASTQ

Agora vamos testar a conversão de arquivos. Primeiro, vamos converter o nosso arquivo FASTQ em arquivos FASTA e QUAL. Em seguida, faremos a conversão inversa utilizando o módulo Bio.SeqIO.QualityIO.

from Bio import SeqIO
from Bio.SeqIO.QualityIO import PairedFastaQualIterator
 
# FASTQ > FASTA
SeqIO.convert("SRR020192.fastq", "fastq", \
"SRR020192.fasta", "fasta")
 
# FASTQ > QUAL
SeqIO.convert("SRR020192.fastq", "fastq", \
"SRR020192.qual", "qual")
 
# FASTQ + QUAL > FASTQ
fastq2 = open("novo_fastq.fastq", "w") 
rec = PairedFastaQualIterator(open("SRR020192.fasta"), \
open("SRR020192.qual"))
i = SeqIO.write(rec, fastq2, "fastq")
fastq2.close()	
print("Foram convertidas %i sequencias FASTA + QUAL em \
formato FASTQ" % i)

.

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