使用字典比较两个 FASTA 文件以获取差异的位置

Question

我有两个文件 f1.fasta 和 f2.fasta。我想比较 f1 和 f2 中的序列，但也想得到核苷酸不同的位置，以便我可以替换它们。

f1 FASTA 示例：

>VFG0127

ATGCCTGGAAATATA...

>VFG0007

TTAGGCATATTTCAT...

f2 FASTA 示例：

>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 不正确，但不确定。祝你好运。

使用字典比较两个 FASTA 文件以获取差异的位置

Comparing two FASTA files to get the position of the difference using dictionary

python

bioinformatics