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

--===============7003055193874941399==
Content-Type: multipart/alternative; boundary="000000000000900505062d8caad5"

--000000000000900505062d8caad5
Content-Type: text/plain; charset="UTF-8"

Good idea, thanks!

I often forget than numpy can work with string based arrays. My initial
working solution (below) was about x10 faster!

Peter


import timeit
import numpy as np


def counts(query_seq, subject_seq):
    gap = "-"
    matches = non_gap_mismatches = either_gapped = 0
    for q, s in zip(query_seq, subject_seq, strict=True):
        # Cache these two as booleans:
        q_gap = q == gap
        s_gap = s == gap
        if q_gap and s_gap:
            continue  # does not contribute to aln_length
        if q == s:
            matches += 1
        elif q_gap or s_gap:
            # Don't want to count this towards coverage (max 100%)
            either_gapped += 1
        else:
            non_gap_mismatches += 1
    return matches, non_gap_mismatches, either_gapped


def np_counts(q_array, s_array):
    q_non_gaps = q_array != b"-"
    s_non_gaps = s_array != b"-"
    naive_matches = q_array == s_array  # includes double gaps!
    matches = int((naive_matches & q_non_gaps).sum())
    one_gapped = q_non_gaps ^ s_non_gaps
    non_gap_mismatches = int((~naive_matches & ~one_gapped).sum())
    either_gapped = int(one_gapped.sum())
    return matches, non_gap_mismatches, either_gapped


SEQ1 = "ACGT-AACCGATG-ATCGTATCGTAG-CGCGTATATGGGG--TTTCGT" * 100000
SEQ2 = "ACGT-AGGTAGCATGCTA-G-CCGTAGCGCGTATATGGGG--TTACGT" * 100000
SEQ1a = np.array(list(SEQ1), "S1")  # bytes!
SEQ2a = np.array(list(SEQ2), "S1")  # bytes!
assert counts(SEQ1, SEQ2) == np_counts(SEQ1a, SEQ2a)
print(timeit.timeit("counts(SEQ1, SEQ2)", number=100, globals=globals()))
print(timeit.timeit("np_counts(SEQ1a, SEQ2a)", number=100,
globals=globals()))

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

<div dir=3D"ltr"><div>Good idea, thanks!</div><div><br></div><div>I often f=
orget than numpy can work with string based arrays. My initial working solu=
tion (below) was about x10 faster!</div><div><br></div><div>Peter</div><div=
><br></div><div><br></div><div>import timeit<br>import numpy as np<br><br><=
br>def counts(query_seq, subject_seq):<br>=C2=A0 =C2=A0 gap =3D &quot;-&quo=
t;<br>=C2=A0 =C2=A0 matches =3D non_gap_mismatches =3D either_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 # Cache these two as booleans:<br>=C2=A0 =C2=
=A0 =C2=A0 =C2=A0 q_gap =3D q =3D=3D gap<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 s_g=
ap =3D s =3D=3D gap<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 if q_gap and s_gap:<br>=
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 continue =C2=A0# does not contrib=
ute to aln_length<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 if 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 elif q_gap or s_gap:<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 # =
Don&#39;t want to count this towards coverage (max 100%)<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 else:<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 non_gap_mismatche=
s +=3D 1<br>=C2=A0 =C2=A0 return matches, non_gap_mismatches, either_gapped=
<br><br><br>def np_counts(q_array, s_array):<br>=C2=A0 =C2=A0 q_non_gaps =
=3D q_array !=3D b&quot;-&quot;<br>=C2=A0 =C2=A0 s_non_gaps =3D s_array !=
=3D b&quot;-&quot;<br>=C2=A0 =C2=A0 naive_matches =3D q_array =3D=3D s_arra=
y =C2=A0# includes double gaps!<br>=C2=A0 =C2=A0 matches =3D int((naive_mat=
ches &amp; q_non_gaps).sum())<br>=C2=A0 =C2=A0 one_gapped =3D q_non_gaps ^ =
s_non_gaps<br>=C2=A0 =C2=A0 non_gap_mismatches =3D int((~naive_matches &amp=
; ~one_gapped).sum())<br>=C2=A0 =C2=A0 either_gapped =3D int(one_gapped.sum=
())<br>=C2=A0 =C2=A0 return matches, non_gap_mismatches, either_gapped<br><=
br><br>SEQ1 =3D &quot;ACGT-AACCGATG-ATCGTATCGTAG-CGCGTATATGGGG--TTTCGT&quot=
; * 100000<br>SEQ2 =3D &quot;ACGT-AGGTAGCATGCTA-G-CCGTAGCGCGTATATGGGG--TTAC=
GT&quot; * 100000<br>SEQ1a =3D np.array(list(SEQ1), &quot;S1&quot;) =C2=A0#=
 bytes!<br>SEQ2a =3D np.array(list(SEQ2), &quot;S1&quot;) =C2=A0# bytes!<br=
>assert counts(SEQ1, SEQ2) =3D=3D np_counts(SEQ1a, SEQ2a)<br>print(timeit.t=
imeit(&quot;counts(SEQ1, SEQ2)&quot;, number=3D100, globals=3Dglobals()))<b=
r>print(timeit.timeit(&quot;np_counts(SEQ1a, SEQ2a)&quot;, number=3D100, gl=
obals=3Dglobals()))<br><br></div></div>

--000000000000900505062d8caad5--

--===============7003055193874941399==
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

--===============7003055193874941399==--