在不匹配("N")的引用上查找字符串(DNA序列)的位置

Question

我试图沿着基因组对齐找到起始和结束位置，因为启动区域是不连续的，所以本质上有两个区域。下面是一个简化的例子：

genome = "GCTAGCTAGCTAGCTGACGCGATCGATGCTAGTTTCGGTGCGCTtAAAAAAGCTTGCATGAAGGCTAGCTA"
primer = "AGCTGANNNNNTCGATGC"

这将像这样对齐：

​

​

我需要找到这个启动区在基因组上的起始和结束位置。

我试着简单地将N拆分为字符串，如下所示：

regions = primer.split('N')
start = genome.find(regions[0])
end = genome.find(regions[-1]) + len(regions[-1])

这样做的问题是，在大的基因组比对中，经常会有较短区域的重复，所以我最终得到了错误的位置。据我所知，BioPython中没有任何东西可以做到这一点，而且pairwise2也没有返回开始和结束位置的方法。

谢谢。

Answer 1

在没有正则表达式的情况下，你可以用这种方式解决这个问题：

genome = "GCTAGCTAGCTAGCTGACGCGATCGATGCTAGTTTCGGTGCGCTtAAAAAAGCTTGCATGAAGGCTAGCTA"
primer = "AGCTGANNNNNTCGATGC"

prefix, suffix = primer[:primer.find("N")], primer[primer.rfind("N") + 1:]

start_pos = end_pos = -1
while True:
    start_pos = genome.find(prefix, start_pos + 1)
    end_pos = start_pos < 0 and start_pos or genome.find(suffix, start_pos + len(prefix)) + len(suffix)
    if start_pos < 0 or end_pos - start_pos == len(primer):
       break

print(start_pos, end_pos)

使用正则表达式：

import re
...

pattern = re.compile(primer.replace("N", "."))
match = pattern.search(genome)
if match:
    start_pos, end_pos = match.span()
    print(start_pos, end_pos)

要以问题中的格式打印，请使用下一步代码：

print(genome)
print(" " * start_pos + "|" * len(prefix) + "." * (len(primer) - len(suffix) - len(prefix)) + "|" * len(suffix))
print("-" * start_pos + primer + "-" * (len(genome) - end_pos))

Answer 2

您也可以简单地使用regex来实现这一点。

import re

genome = "GCTAGCTAGCTAGCTGACGCGATCGATGCTAGTTTCGGTGCGCTtAAAAAAGCTTGCATGAAGGCTAGCTA"
primer = "AGCTGANNNNNTCGATGC"

#create regular expression from primer
reg = '[A-Z]+'.join([i for i in primer.split('N') if len(i)>0])

#Search for index where regular expression matches
idx = re.search(reg, genome).span()

#Output = (11,29) which is start and end of the match

#print both with primer aligned based on start to end index of match
print(genome)
print(''.join((['_']*idx[0])+list(primer)+(['_']*(len(genome)-idx[1]))))

GCTAGCTAGCTAGCTGACGCGATCGATGCTAGTTTCGGTGCGCTtAAAAAAGCTTGCATGAAGGCTAGCTA
___________AGCTGANNNNNTCGATGC__________________________________________

在regex101.com上测试的正则表达式

​

​

Answer 3

import re

genome = "GCTAGCTAGCTAGCTGACGCGATCGATGCTAGTTTCGGTGCGCTtAAAAAAGCTTGCATGAAGGCTAGCTA"
primer = "AGCTGANNNNNTCGATGC"

split = re.split('N+', primer)
middle = len(primer)-len(''.join(split))

r = f'{split[0]}\w{{{middle}}}{split[-1]}'

re.finditer(r,genome)

# Taking index 0 here assuming there will be one match
# If there are multiple matches handle the printing by looping perhaps
pos = [(m.start(0), m.end(0)) for m in re.finditer(r, genome)][0]

print(genome)
print(' '*pos[0] + '|'*len(split[0]) + '.'*middle + '|'*len(split[-1]))
print('-'*pos[0] + primer + '-'*int(len(genome)-int(pos[0]+len(primer))))

输出

GCTAGCTAGCTAGCTGACGCGATCGATGCTAGTTTCGGTGCGCTtAAAAAAGCTTGCATGAAGGCTAGCTA
           ||||||.....|||||||
-----------AGCTGANNNNNTCGATGC------------------------------------------