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

Question

我试图获得一系列序列中包含的插入和间隔的数量，这些序列与它们被对齐的参考有关；因此，所有序列现在都具有相同的长度。

例如

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

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

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

^{pr2}$

我有点像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)

输出

^{pr2}$

---

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

^{3}$

输出

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}

Answer 2

使用生物疗法和裸体：

from Bio import AlignIO
from collections import Counter
import numpy as np


alignment = AlignIO.read("alignment.fasta", "fasta")

events = []

for i in range(alignment.get_alignment_length()):
    this_column = alignment[:, i]

    # Mark insertions, polymorphism and deletions following PM 2Ring notation
    events.append(["U" if b == this_column[0] else
                   "I" if this_column[0] == "-" else
                   "P" if b != "-" else
                   "D" for b in this_column])

# Apply a Counter over the columns (axis 0) of the array
print(np.apply_along_axis(Counter, 0, np.array(events)))

这应输出与对齐顺序相同的计数数组：

^{pr2}$