我想通过一个基因,从每个feature.type =='mRNA'获得一个包含外显子/内含子边界的10 go长序列的列表。我似乎需要使用compoundLocation,以及'join'中使用的位置,但我不知道如何做到这一点,也无法找到教程。
有谁能给我举个例子或者给我介绍一个教程吗?
发布于 2015-01-04 18:07:21
假设所有的信息都以您在评论中所显示的确切格式显示,并且您希望在每个引论/外显子边界的任何一侧查找20 bp,这样的情况可能是一个开始:
编辑:--如果您实际上是从GenBank记录开始的,那么就不难了。假设您要查找的完整连接字符串位于CDS功能信息中,那么:
for f in record.features:
if f.type == 'CDS':
jct_info = str(f.location)将“位置”信息转换为字符串,您可以继续如下所示。
(有一些方法可以直接处理位置信息,而无需转换为字符串--特别是可以使用“提取”将拼接序列直接从父序列中提取出来--但是,通过转换为str然后int,您想要做的步骤更快、更容易完成。)
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]']现在,您可以循环遍历start:end值的列表,将它们从文本中提取出来,并将它们转换为int,这样就可以将它们用作序列索引。就像这样:
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)https://stackoverflow.com/questions/27765857
复制相似问题