比較兩個FASTA文件以使用字典獲得差異的位置

Question

我有兩個文件f1.fasta和f2.fasta。我想比較f1和f2中的序列，但也可以得到核苷酸不同的位置，以便我可以替換它們。

例 f1 FASTA

：

例 f2 FASTA

>VFG0127 

ATGCCTGGAAATATA... 

>VFG0007 

TTAGGCATATTTCAT... 

：

>VFG0127 

ATGCCTGGXXXTATA... 

>VFG0007 

TTAJGCATATSTCAT... 

我想獲得，例如：VFG0127 |位置7，X應該是A ...

我試過這段代碼，但我沒有去任何地方

dict_1 = {} 
dict_2 = {} 

with open(f1, 'r') as f1, open (f2, 'r') as f2: 
    for line in f1: 
     if line.startswith('>'): 
      id_acc1 = line.strip() 
      seq_1 = f1.next().strip() 
      dict_1[id_acc1]=seq_1 
      #print dict_1 
    for line in f2: 
     if line.startswith('>'): 
      id_acc2 = line.strip() 
      seq_2 = f2.next().strip() 
      dict_2[id_acc2]=seq_2 
      #print dict_2 

    diffkeys = [k for k in dict_1.values()[index] if dict_1[k] != dict_2[k]] 
    for k in diffkeys: 
     print k, ':', dict_1[k], '->', dict_2[k] 

我在這個問題上花了數小時，我無法使它工作。 請我仍然是一個初學者，一個簡單的代碼將不勝感激。

Answer 1

我認爲這應該工作。將取決於你的fasta文件的確切程度。

f1 = open("f1.fasta","r").readlines() 
f2 = open("f2.fasta","r").readlines() 

## Read the files 
dict1 = {} 
dict2 = {} 
currentID = "" 
for l in f1: 
    line = l.strip() 
    if line[0] == ">": 
     currentID = line[1:] 
     dict1[line[1:]] = "" 
    else: 
     dict1[currentID] = dict1[currentID]+line 

currentID = "" 

for l in f2: 
    line = l.strip() 
    if line[0] == ">": 
     currentID = line[1:] 
     dict2[line[1:]] = "" 
    else: 
     dict2[currentID] = dict2[currentID]+line 

##Assuming that both sequences have same length 
for key in dict1.keys(): 
    if dict1[key] != dict2[key]: 
     for i in range(len(dict1[key])): 
      if dict1[key][i] != dict2[key][i]: 
       print(key, i, dict1[key][i], dict2[key][i]) 

我認爲你的diffkeys是不正確的，但不知道。好運。