How to get the number of non-shared insertions and gaps in a sequence in Python?
<p>I am trying to obtain the number of insertions and gaps contained in a series of sequences with relation to a reference with which they were aligned; therefore, <strong>all sequences are now of the same length.</strong></p>

<p>For instance</p>

<pre><code>&gt;reference
AGCAGGCAAGGCAA--GGAA-CCA
&gt;sequence1
AAAA---AAAGCAATTGGAA-CCA
&gt;sequence2
AGCAGGCAAAACAA--GGAAACCA
</code></pre>

<p>In this example, sequence1 has two insertions (two T) and three gaps. The last gap should not be counted since it appears both in the reference and sequence1. Sequence2 has one insertion (an A before the last triplet) and  no gaps. (Again, the gaps are shared with the reference and should no enter in the count.). There are also 3 polymorphisms in sequence 1 and 2 in sequence 2.</p>

<p>My current script is able to give an estimate of the differences but not the count of "relevant gaps and insertions" as described above. For example</p>

<pre><code>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
</code></pre>

<p>I am kind of a Python newy and still learning the tools of this language. Is there a way to achieve this? I am thinking about spliting the sequences and iteratively compare one position at a time and count the difference but I am not sure if it possible in Python (not to mention that it would be horribly slow.)</p>

<p>Using Biopython and numpy:</p>

<pre><code>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)))
</code></pre>

<p>This should output an array of Counts in the same order as the alignment:</p>

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

python bioinformatics biopython fasta