python - 如何在Python中获取序列中非共享插入和间隙的数量？

Question

我正在尝试获取一系列序列中包含的插入和间隙的数量，以及与它们对齐的引用的关系；因此，所有序列现在都具有相同的长度。

例如

>reference
AGCAGGCAAGGCAA--GGAA-CCA
>sequence1
AAAA---AAAGCAATTGGAA-CCA
>sequence2
AGCAGGCAAAACAA--GGAAACCA

在此示例中，sequence1 有两次插入(两个 T)和三个间隙。最后一个间隙不应计算在内，因为它同时出现在引用和序列 1 中。 Sequence2 有一个插入(最后一个三元组之前有一个 A)并且没有间隙。 (同样，间隙与引用共享，不应输入计数中。)。序列1中也有3个多态性，序列2中也有2个多态性。

我当前的脚本能够给出差异的估计，但不能给出如上所述的“相关间隙和插入”的计数。例如

records = list(SeqIO.parse(file("sequences.fasta"),"fasta"))
reference = records[0] #reference is the first sequence in the file
del records[0]

for record in records:
   gaps = record.seq.count("-") - reference.seq.count("-")
   basesinreference = reference.seq.count("A") + reference.seq.count("C") + reference.seq.count("G") + reference.seq.count("T")
   basesinsequence = record.seq.count("A") + record.seq.count("C") + record.seq.count("G") + record.seq.count("T")
   print(record.id)
   print(gaps)
   print(basesinsequence - basesinreference)
#Gives
sequence1
1 #Which means sequence 1 has one more Gap than the reference
-1 #Which means sequence 1 has one base less than the reference
sequence2
-1 #Which means sequence 2 has one Gap less than the reference
1 #Which means sequence 2 has one more base than the reference

我是一个 Python 新手，仍在学习这种语言的工具。有办法实现这一点吗？我正在考虑分割序列并一次迭代地比较一个位置并计算差异，但我不确定在 Python 中是否可能(更不用说它会非常慢。)

Answer 1

这是 zip 函数的工作。我们并行迭代引用序列和测试序列，查看其中一个在当前位置是否包含 -。我们使用该测试的结果来更新字典中插入、删除和未更改的计数。

def kind(u, v):
    if u == '-':
        if v != '-':
            return 'I'  # insertion
    else:
        if v == '-':
            return 'D'  # deletion
    return 'U'          # unchanged

reference = 'AGCAGGCAAGGCAA--GGAA-CCA'

sequences = [
    'AGCA---AAGGCAATTGGAA-CCA',
    'AGCAGGCAAGGCAA--GGAAACCA',
]

print('Reference')
print(reference)
for seq in sequences:
    print(seq)
    counts = dict.fromkeys('DIU', 0)
    for u, v in zip(reference, seq):
        counts[kind(u, v)] += 1
    print(counts)

输出

Reference
AGCAGGCAAGGCAA--GGAA-CCA
AGCA---AAGGCAATTGGAA-CCA
{'I': 2, 'D': 3, 'U': 19}
AGCAGGCAAGGCAA--GGAAACCA
{'I': 1, 'D': 0, 'U': 23}

---

这是一个更新版本，它还检查多态性。

def kind(u, v):
    if u == '-':
        if v != '-':
            return 'I'  # insertion
    else:
        if v == '-':
            return 'D'  # deletion
        elif v != u:
            return 'P'  # polymorphism
    return 'U'          # unchanged

reference = 'AGCAGGCAAGGCAA--GGAA-CCA'

sequences = [
    'AAAA---AAAGCAATTGGAA-CCA',
    'AGCAGGCAAAACAA--GGAAACCA',
]

print('Reference')
print(reference)
for seq in sequences:
    print(seq)
    counts = dict.fromkeys('DIPU', 0)
    for u, v in zip(reference, seq):
        counts[kind(u, v)] += 1
    print(counts)

输出

Reference
AGCAGGCAAGGCAA--GGAA-CCA
AAAA---AAAGCAATTGGAA-CCA
{'D': 3, 'P': 3, 'I': 2, 'U': 16}
AGCAGGCAAAACAA--GGAAACCA
{'D': 0, 'P': 2, 'I': 1, 'U': 21}