Допустим, у меня есть идентификатор GenBank для генома (BA000007.3) и список с координатами генома для блоков:
blocks = [(1443423, 1444553), ...]
Мне нужно выяснить, какие гены (или частикакие гены ) покрыты этими блоками.
Я пробовал это:
from Bio import Entrez, SeqIO
Entrez.email = "some@example.com"
with Entrez.efetch(db="nucleotide", rettype="gb", retmode="text", id="BA000007.3") as handle:
seq_record = SeqIO.read(handle, "gb")
sub_record = seq_record[1443423:1444553]
print(len(sub_record.features))
Возвращает 0, поскольку гены полностью не покрыты этим блоком.Однако, если вы откроете запись NCBI в браузере и измените показанный регион, вы увидите, что есть 2 гена, которые частично покрываются этим блоком.
Я также пытался вычленять с помощью SeqFeature и FeatureLocation, но он тоже не работал (новая последовательность вообще не имеет функций).
Так есть ли способ сделать этобез ручной проверки каждого блока для каждого генома?