在 fasta 序列中查找核苷酸子序列答案

【问题标题】：Find nucleotide subsequence within fasta sequences在 fasta 序列中查找核苷酸子序列
【发布时间】：2016-11-30 00:27:16
【问题描述】：

我必须编写一个函数，输入一个包含带有ambiguous symbols (IUPAC) 的 dna 序列的 FASTA 文件。给定 FASTA 文件的名称和明确的 DNA 字符串，我想写出给定序列可能是子序列的序列标识符（'>' 标题）。我想在不生成所有可能的序列的情况下实现这一点，并且子序列可以有不明确的符号以及 FASTA 文件中的序列。示例：序列“ACC”可以是“CGMBHTW”的子序列。有人可以帮我吗？

【问题讨论】：

你的尝试是什么？ -> 到目前为止显示您的代码
您能否提供任何代码让我们开始使用？至于我，我不知道这个问题背后的生物学，所以更清楚你要解决的问题也会有所帮助..
您能否提供一些测试示例来显示每个测试的输入和正确输出？
Can Biopython perform Seq.find() accounting for ambiguity codes的可能重复

标签： python bioinformatics fasta

【解决方案1】：

您可以尝试将“通用”核苷酸定义为它所代表的一组字母，然后将您的序列转换为此类集合的列表，并扫描一个序列是否存在与另一个序列兼容的位置。

这可能不是最有效的代码，但似乎可以工作（请仔细检查我的循环索引是否正确......）。

A = {"A"}
C = {"C"}
G = {"G"}
T = {"T"}
R = A | G
Y = C | T
S = G | C
W = A | T
K = G | T
M = A | C
B = C | G | T
D = A | G | T
H = A | C | T
V = A | C | G
N = {"A", "C", "G", "T"}

letter2nucl = {
    "A" : A,
    "C" : C,
    "G" : G,
    "T" : T,
    "R" : R,
    "Y" : Y,
    "S" : S,
    "W" : W,
    "K" : K,
    "M" : M,
    "B" : B,
    "D" : D,
    "H" : H,
    "V" : V,
    "N" : N}

def is_subseq(seq1, seq2):
    l1 = len(seq1)
    l2 = len(seq2)
    nucls1 = [letter2nucl[letter] for letter in seq1]
    nucls2 = [letter2nucl[letter] for letter in seq2]
    i = 0
    while i < 1 + l2 - l1:
        subseq = True
        for j, nucl in enumerate(nucls1):
            if not (nucls2[i+j] & nucl):
                # empty set intersection
                subseq = False
                break
        if subseq:
            return True
        i += 1
    return False

seq1 = "ACC"
seq2 = "CGMBHTW"

if is_subseq(seq1, seq2):
    print("%s is subsequence of %s" % (seq1, seq2))

seq1 = "GRT"
seq2 = "AATCBAT"

if is_subseq(seq1, seq2):
    print("%s is subsequence of %s" % (seq1, seq2))

结果是：

ACC is subsequence of CGMBHTW
GRT is subsequence of AATCBAT

然后，您可以在使用 Biopython 的 SeqIO 函数读取的序列上使用它。

【讨论】：