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>This is a job for the <code>zip</code> function. We iterate over the reference and a test sequence in parallel, seeing if either one contains a <code>-</code> at the current position. We use the result of that test to update counts of insertions, deletions and unchanged in a dictionary.</p>

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

<p><strong>output</strong></p>

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

<hr>

<p>Here's an updated version that also checks for polymorphism.</p>

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

<p><strong>output</strong></p>

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

python bioinformatics biopython fasta