<div dir="ltr"><div>Good idea, thanks!</div><div><br></div><div>I often forget than numpy can work with string based arrays. My initial working solution (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>    gap = "-"<br>    matches = non_gap_mismatches = either_gapped = 0<br>    for q, s in zip(query_seq, subject_seq, strict=True):<br>        # Cache these two as booleans:<br>        q_gap = q == gap<br>        s_gap = s == gap<br>        if q_gap and s_gap:<br>            continue  # does not contribute to aln_length<br>        if q == s:<br>            matches += 1<br>        elif q_gap or s_gap:<br>            # Don't want to count this towards coverage (max 100%)<br>            either_gapped += 1<br>        else:<br>            non_gap_mismatches += 1<br>    return matches, non_gap_mismatches, either_gapped<br><br><br>def np_counts(q_array, s_array):<br>    q_non_gaps = q_array != b"-"<br>    s_non_gaps = s_array != b"-"<br>    naive_matches = q_array == s_array  # includes double gaps!<br>    matches = int((naive_matches & q_non_gaps).sum())<br>    one_gapped = q_non_gaps ^ s_non_gaps<br>    non_gap_mismatches = int((~naive_matches & ~one_gapped).sum())<br>    either_gapped = int(one_gapped.sum())<br>    return matches, non_gap_mismatches, either_gapped<br><br><br>SEQ1 = "ACGT-AACCGATG-ATCGTATCGTAG-CGCGTATATGGGG--TTTCGT" * 100000<br>SEQ2 = "ACGT-AGGTAGCATGCTA-G-CCGTAGCGCGTATATGGGG--TTACGT" * 100000<br>SEQ1a = np.array(list(SEQ1), "S1")  # bytes!<br>SEQ2a = np.array(list(SEQ2), "S1")  # bytes!<br>assert counts(SEQ1, SEQ2) == np_counts(SEQ1a, SEQ2a)<br>print(timeit.timeit("counts(SEQ1, SEQ2)", number=100, globals=globals()))<br>print(timeit.timeit("np_counts(SEQ1a, SEQ2a)", number=100, globals=globals()))<br><br></div></div>