Как получить количество несвязанных вставок и пробелов в последовательности в Python?

Question

Я пытаюсь получить количество вставок и зазоров, содержащихся в последовательности последовательностей, относительно ссылки, с которой они были выровнены; поэтому все последовательности имеют одинаковую длину.

Например

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

В этом примере sequence1 имеет две вставки (два Т) и трех пробелов. Последний пробел не следует учитывать, поскольку он появляется как в ссылке, так и в последовательности1. Последовательность2 имеет одну вставку (A перед последним триплетом) и отсутствие пробелов. (Опять же, пробелы разделяются ссылкой и не должны вводить счет.). Существует также 3 полиморфизма в последовательностях 1 и 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 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})]]