[Biopython] Comparing two already aligned sequences quickly to count gaps/matches/mismatches

--===============6076023287401778408==
Content-Type: multipart/alternative; boundary="000000000000aaa1fc062ceb9640"

--000000000000aaa1fc062ceb9640
Content-Type: text/plain; charset="UTF-8"

Hello all, and Michiel in particular,

I am wondering if any of the pairwise alignment code in Bio.Align (much of
which is written in C for speed) could help with this use case?:

I have lots of pairs of pre-aligned sequences (imported from an external
MSA file), for which I am doing something like this:

```python
def count_matches_etc(query_seq, subject_seq):
    assert len(query_seq) == len(subject_seq), "Should be same length"
    matches = non_gap_mismatches = either_gapped = both_gapped = 0
    for q, s in zip(query_seq, subject_seq, strict=True):
        if q == "-" and s == "-":
            both_gapped += 1
        elif q == "-" or s == "-":
            either_gapped += 1
        elif q == s:
            matches += 1
        else:
            non_gap_mismatches += 1
    assert matches + non_gap_mismatches + either_gapped + both_gapped ==
len(query_seq)
    return matches, non_gap_mismatches, either_gapped, both_gapped


# Test case
assert (9, 1, 2, 1) == count_matches_etc("ACGTAC-TAC-GT", "AGGT-CGTAC-GT")
```

Sticking with Python that could be optimized (e.g. I am currently using
this with sequences of a million base pairs but few gaps), however I have
written this example with clarity foremost in mind.

Thank you,

Peter

--000000000000aaa1fc062ceb9640
Content-Type: text/html; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable

<div dir=3D"ltr">Hello all, and Michiel in particular,<br><br>I am wonderin=
g if any of the pairwise alignment code in Bio.Align (much of which is writ=
ten in C for speed) could help with this use case?:<br><br>I have lots of p=
airs of pre-aligned sequences (imported from an external MSA file), for whi=
ch I am doing something like this:<br><br>```python<br>def count_matches_et=
c(query_seq, subject_seq):<br>=C2=A0 =C2=A0 assert len(query_seq) =3D=3D le=
n(subject_seq), &quot;Should be same length&quot;<br>=C2=A0 =C2=A0 matches =
=3D non_gap_mismatches =3D either_gapped =3D both_gapped =3D 0<br>=C2=A0 =
=C2=A0 for q, s in zip(query_seq, subject_seq, strict=3DTrue):<br>=C2=A0 =
=C2=A0 =C2=A0 =C2=A0 if q =3D=3D &quot;-&quot; and s =3D=3D &quot;-&quot;:<=
br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 both_gapped +=3D 1<br>=C2=A0 =
=C2=A0 =C2=A0 =C2=A0 elif q =3D=3D &quot;-&quot; or s =3D=3D &quot;-&quot;:=
<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 either_gapped +=3D 1<br>=C2=
=A0 =C2=A0 =C2=A0 =C2=A0 elif q =3D=3D s:<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =
=C2=A0 =C2=A0 matches +=3D 1<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 else:<br>=C2=A0=
 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 non_gap_mismatches +=3D 1<br>=C2=A0 =C2=
=A0 assert matches + non_gap_mismatches + either_gapped + both_gapped =3D=
=3D len(query_seq)<br>=C2=A0 =C2=A0 return matches, non_gap_mismatches, eit=
her_gapped, both_gapped<br><br><br># Test case<br>assert (9, 1, 2, 1) =3D=
=3D count_matches_etc(&quot;ACGTAC-TAC-GT&quot;, &quot;AGGT-CGTAC-GT&quot;)=
<br>```<br><br>Sticking with Python that could be optimized (e.g. I am curr=
ently using this with sequences of a million base pairs but few gaps), howe=
ver I have written this example with clarity foremost in mind.<br><div><br>=
</div><div>Thank you,</div><div><br></div>Peter</div>

--000000000000aaa1fc062ceb9640--

--===============6076023287401778408==
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

_______________________________________________
Biopython mailing list  -  [email protected]
https://mailman.open-bio.org/mailman/listinfo/biopython

--===============6076023287401778408==--