es.davy.ai

Preguntas y respuestas de programación confiables

¿Tienes una pregunta?

Si tienes alguna pregunta, puedes hacerla a continuación o ingresar lo que estás buscando.

Biopython | ¿Por qué solo obtengo valores de “e” de 0.0?

12 Jul, 2023 Programación

Estoy buscando obtener los 5 mejores resultados de alineación de blast para la secuencia de ARNr 16S que obtuve con el módulo Entrez. Sin embargo, al analizar mis alineaciones, todos los puntajes de valor e (e-value) que obtengo de record.alignment.hsp.expect son exactamente 0.0.

Solo con inspeccionar visualmente las alineaciones, las secuencias no eran 100% iguales.

from Bio import SeqIO
from Bio.Seq import Seq
from Bio import Entrez
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

Entrez.email = "obsto123@gmail.com"
handle = Entrez.efetch(db="nucleotide", id="NC_000913.1", rettype="gb", retmode="text")
seq_record = SeqIO.read(handle, "genbank")

records = []
for feature in seq_record.features:
    if feature.type == "rRNA" and "16S" in feature.qualifiers["product"][0]:
        records.append(seq_record.seq[feature.location.start.position:feature.location.end.position])
rRNA = records[0]

print("Ahora realizando búsqueda de blast...")

blast_data = NCBIWWW.qblast("blastn", "nt", records[0])

with open("blast_xml.xml", "w") as f:
    f.write(blast_data.read())
xml = open("blast_xml.xml", "r")
blast_record = NCBIXML.read(xml)

for alignment in blast_record.alignments:
    for hsp in alignment.hsps:
        print(hsp.expect)
Tags:  , , ,

Answer

  1. Avatar for davy.ai
    davy.ai

    Solución

    El problema de obtener todos los puntajes de e-value como 0.0 se debe a utilizar la variable incorrecta para obtener la secuencia para realizar el blasting. En lugar de utilizar records[0], que es solo la secuencia del 16S rRNA, se debe utilizar la secuencia completa seq_record.seq de la siguiente manera:

    blast_data = NCBIWWW.qblast("blastn", "nt", seq_record.seq)

    Con este cambio, los puntajes de e-value deberían calcularse correctamente para cada alineación. Para obtener las 5 mejores alineaciones basadas en el puntaje de e-value, se puede utilizar la función sorted incorporada de Python con el parámetro key establecido en el valor “expect” de cada alineación de la siguiente manera:

    top_results = sorted(blast_record.alignments, key=lambda x: x.hsps[0].expect)[:5]

    Aquí, se utiliza la función lambda para extraer el valor “expect” de cada alineación y ordenar los resultados en orden ascendente. El [:5] al final se utiliza para seleccionar las primeras 5 alineaciones con los valores de e-value más bajos. El código final podría verse así:

    from Bio import SeqIO
from Bio import Entrez
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

Entrez.email = "obsto123@gmail.com"
handle = Entrez.efetch(db="nucleotide", id="NC_000913.1", rettype="gb", retmode="text")
seq_record = SeqIO.read(handle, "genbank")

print("Ahora realizando la búsqueda de blasting...")

blast_data = NCBIWWW.qblast("blastn", "nt", seq_record.seq)

with open("blast_xml.xml", "w") as f:
    f.write(blast_data.read())

xml = open("blast_xml.xml", "r")
blast_record = NCBIXML.read(xml)

top_results = sorted(blast_record.alignments, key=lambda x: x.hsps[0].expect)[:5]

for i, alignment in enumerate(top_results):
    print(f"\nAlineación {i+1}:\nSujeto: {alignment.title}")
    print(f"Puntaje de e-value: {alignment.hsps[0].expect}")
    print(f"Longitud: {alignment.length}")
    print(f"Puntaje: {alignment.hsps[0].score}")
    print(f"Huecos: {alignment.hsps[0].gaps}")
    print(f"Consulta: {alignment.hsps[0].query}")
    print(f"Coincidencia: {alignment.hsps[0].match}")
    print(f"Sujeto: {alignment.hsps[0].sbjct}")

    Esto debería imprimir los detalles de las 5 mejores alineaciones con los valores de e-value más bajos, incluido el título de la alineación, e-value, longitud, puntaje, huecos y las alineaciones de secuencia de ADN para la consulta, coincidencia y sujeto.

Comments are closed.