我想遍历一个基因并获得一个 10bp 长的序列列表,其中包含来自每个 feature.type =='mRNA'
的外显子/内含子边界。似乎我需要使用 compoundLocation
,以及 'join'
中使用的位置,但我不知道该怎么做,也找不到教程。
谁能给我举个例子或指点教程吗?
最佳答案
假设所有信息都以您在评论中显示的确切格式显示,并且您正在寻找每个介绍/外显子边界两侧的 20 bp,这样的事情可能是一个开始:
编辑: 如果您实际上是从 GenBank 记录开始的,那么这并不会更难。假设您要查找的完整联结字符串在 CDS 功能信息中,则:
for f in record.features:
if f.type == 'CDS':
jct_info = str(f.location)
将“位置”信息转换为字符串,您可以继续执行以下操作。
(有一些方法可以直接处理位置信息而无需转换为字符串 - 特别是您可以使用“提取”将拼接序列直接从父序列中拉出 - 但是您想要的步骤通过先转换为 str 再转换为 int,do 会更快更容易地完成。)
import re
jct_info = "join{[0:229](+), [11680:11768](+), [11871:12135](+), [15277:15339](+), [16136:16416](+), [17220:17471](+), [17547:17671](+)"
jctP = re.compile("\[\d+\:\d+\]")
jcts = jctP.findall(jct_info)
jcts
['[0:229]', '[11680:11768]', '[11871:12135]', '[15277:15339]', '[16136:16416]', '[17220:17471]', '[17547:17671]']
现在您可以循环遍历开始:结束值列表,将它们从文本中提取出来并将它们转换为整数,以便您可以将它们用作序列索引。像这样:
for jct in jcts:
(start,end) = jct.replace('[', '').replace(']', '').split(':')
try: # You need to account for going out of index, e.g. where start = 0
start_20_20 = seq[int(start)-20:int(start)+20]
except IndexError:
# do your alternatives e.g. start = int(start)
关于python - 寻找基因中的外显子/内含子边界,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27765857/