fasta配列内のヌクレオチド部分配列を見つける

Question

ambiguous symbols (IUPAC)のdna配列を含むFASTAファイルを入力して関数を書く必要があります。 FASTAファイルの名前とあいまいなDNA文字列が与えられているので、与えられたシーケンスがサブシーケンスである可能性のあるシーケンス（'>'ヘッダ）の識別子を書きたいと思います。私は可能なすべてのシーケンスを生成せずにこれを実装したいと思います。サブシーケンスはFASTAファイル内のシーケンスと同様にあいまいなシンボルを持つことがあります。例：シーケンス "ACC"は "CGMBHTW"のサブシーケンスである可能性があります。 誰かが私を助けることができますか？

Answer 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機能を使用して読み取ります。