查找fasta序列內的核苷酸子序列

我必須編寫一個函數，輸入一個包含dna序列的FASTA文件，其中ambiguous symbols (IUPAC)。鑑於FASTA文件的名稱和明確的DNA字符串，我想寫出序列的標識符（'>'標題），其中給定的序列可能是子序列。我希望在不生成所有可能的序列的情況下執行此操作，並且子序列可能具有模糊符號以及FASTA文件中的序列。例如：序列「ACC」可能是「CGMBHTW」的子序列。有人可以幫我嗎？查找fasta序列內的核苷酸子序列

來源

2016-11-30 H.F.S C.

你的嘗試是什麼？ - >顯示你的代碼到目前爲止 – shash678

你能提供*任何*代碼給我們開始使用嗎？至於我，我不知道這個問題背後的生物學，所以你想要解決的問題的更清晰也將有所幫助。 – davedwards

你能否提供一些測試示例，顯示每個測試的輸入和正確的輸出？ – Steve

您可以嘗試將「廣義」核苷酸定義爲它所表示的一組字母，然後將這些序列轉換爲這些集合的列表，然後掃描一個序列以找到與另一個序列兼容的位置。

下面是可能不是最有效的代碼，但似乎工作（仔細檢查我的循環指數是正確的，雖然...）。

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功能。

來源

2016-12-07 11:05:21 bli

查找fasta序列內的核苷酸子序列

回答

相關問題