如何获取Python中序列中非共享插入和空位的数量？

Question

我试图获得包含在一系列序列中的插入和空位的数量与它们对齐的引用的关系;因此，所有序列现在具有相同的长度。

例如

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

在这个例子中，序列1具有两个插入（2 T）和三个间隙。由于它出现在参考文献和序列1中，因此不应该计算最后的差距。序列2有一个插入（最后一个三联体之前的A）并且没有空位。 （再次，差距与参考共享，不应输入计数。）。在序列1的序列1和2中也有3个多态性。

我的当前脚本能够估计差异，但不能计算上述“相关差距和插入”的数量。例如

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 newy，仍然在学习这种语言的工具。有没有办法做到这一点？我正在考虑拆分序列并迭代比较一个位置并计算差异，但我不确定是否可以在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} 

Answer 2

使用Biopython和numpy的：

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))) 

这应该输出计数的阵列中的相同顺序的取向：

[[Counter({'U': 23}) 
    Counter({'U': 15, 'P': 3, 'D': 3, 'I': 2}) 
    Counter({'U': 21, 'P': 2, 'I': 1})]]