检查字符串中是否存在子序列并允许1-3个不匹配的快速python方法



我想知道是否有一种快速的Python方法可以检查字符串中是否存在子序列,同时允许1-3个不匹配。

例如字符串:";ATGCTGCTGA";子序列";ATGCC";将是可接受的,并返回true。请注意,粗体中有1个不匹配。

我试过使用Bio软件包中的pairwise2功能,但速度很慢。我有1000个字符串,对于每个字符串,我想测试10000个子序列。

这里将优先考虑计算速度。

**注意,我指的不是间隙,而是一个核苷酸(字母A、t、C、G(被另一个取代。

可以压缩字符串,比较左右两边的元组,并计算false。这里false只是1。不应该太慢。。。

st = "ATGCTGCTGA"
s = "ATGCC"
[ x==y for (x,y) in zip(st,s)].count(False)
1

尝试:

ALLOWED_MISMATCHES = 3
s = "ATGCTGCTGA"
subsequence = "ATGCC"
for i in range(len(s) - len(subsequence) + 1):
if sum(a != b for a, b in zip(s[i:], subsequence)) <= ALLOWED_MISMATCHES:
print("Match")
break
else:
print("No Match")

打印:

Match

如果您使用PyPi的正则表达式模块来代替模糊匹配,则正则表达式可以实现这一点:

ATGCC{i<=3:[ATGC]}
  • ATGCC-精确查找"ATGCC">
  • {i<=3:[ATCG]}-允许在核苷酸字符[ATGC]的字符类别内最多插入三个

例如:

import regex as re
s = 'ATGCTGCTGA'
print(bool(re.search(r'ATGCC{i<=3:[ATGC]}', s)))

打印:

True

特别是,如果您的子序列都有相同的长度,您可以尝试k-mer与Biosite匹配,我是该软件包的开发人员。为了允许不匹配,可以使用SimilarityRule在匹配过程中生成类似的子序列。在本例中,设置了规则和替换矩阵,以便枚举具有多达MAX_MISMATCH不匹配的所有子序列。由于k-mer匹配是在C中实现的,它应该运行得很快。

import numpy as np
import biotite.sequence as seq
import biotite.sequence.align as align
K = 5
MAX_MISMATCH = 1

database_sequences = [
seq.NucleotideSequence(seq_str) for seq_str in [
"ATGCTGCTGA",
# ...
]
]
query_sequences = [
seq.NucleotideSequence(seq_str) for seq_str in [
"ATGCC",
# ...
]
]

kmer_table = align.KmerTable.from_sequences(K, database_sequences)
# The alphabet for the substitution matrix
# should not contain unambiguous symbols, such as 'N'
alphabet = seq.NucleotideSequence.alphabet_unamb
matrix_entries = np.full((len(alphabet), len(alphabet)), -1)
np.fill_diagonal(matrix_entries, 0)
matrix = align.SubstitutionMatrix(alphabet, alphabet, matrix_entries)
print(matrix)
print()
# The similarity rule will allow up to MAX_MISMATCH mismatches
similarity_rule = align.ScoreThresholdRule(matrix, -MAX_MISMATCH)
for i, sequence in enumerate(query_sequences):
matches = kmer_table.match(sequence, similarity_rule)
# Colums:
# 0. Sequence position of match in query (first position of k-mer)
# 1. Index of DB sequence in list
# 2. Sequence position of match in DB sequence (first position of k-mer)
print(matches)
print()
index_of_matches = np.unique(matches[:, 1])
for j in index_of_matches:
print(f"Match of subsequence {i} to sequence {j}")
print()

输出:

A   C   G   T
A   0  -1  -1  -1
C  -1   0  -1  -1
G  -1  -1   0  -1
T  -1  -1  -1   0
[[0 0 0]]
Match of subsequence 0 to sequence 0

如果子序列的长度不同,但都很短(长度为10(,则可以为每个子序列长度(K=长度(创建一个KmerTable,并将每个表与具有相应长度的子序列相匹配。如果您的子序列比这个大得多,那么由于KmerTable的内存需求,这种方法可能无法工作。

相关内容

最新更新