How can one concatenate two very large strings and return indices for matches and mismatches?

I have a set of text files that contain two very large character sets of the same length. Symbol sets are DNA sequences, so I will call them seq_1

and seq_2

, and together they are called alignment

. The files look like this:

>HEADER1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>HEADER2
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

      

        
        
        
      

    

Possible characters in sequence 1 below >HEADER1

are ACGTN-

, and possible characters in sequence 2 below >HEADER2

are equal ACGTN-*

.

I want to parse sequences and return two lists of indices, which I will call valid

and mismatch

.

valid

contains all (1-based) indices where the positions in both sequences ("alignment") are in the set ACGT

; mismatch

contains all indices (1 based) where the positions in the alignment are in the set ACGT

but not the same. Thus, mismatch

is a subset of valid

.

The end-point is that I DO NOT increment the index counter at the positions where sequence 1 is "-"

, because they are considered "spaces" in the coordinate system being used.

This example shows the alignment of my expected result:

seq_1      = 'CT-A-CG'  # total length = 5 because of the two gaps
seq_2      = 'CNCTA*G'
valid      = [1,3,5]    # all indices where both sequences are in 'ACGT' without counting gaps
mismatch   = [3]        # subset of valid where alignment does not match

      

        
        
        
      

    

I want to improve on my current code (see below) which involves fetching sequences regularly, encrypting and enumerating non-market sites into a generator object, and then the main step of the time it takes to loop through that generator and populate the two lists. I feel like there must be an array-based solution or solution to this problem itertools

that would be much more efficient than a for loop through a loop and their index zipped together and I'm looking for suggestions.

code:

def seq_divergence(infile):

    re_1 = re.compile('>HEADER1\n([^>]+)', re.MULTILINE)
    re_2 = re.compile('>HEADER2\n([^>]+)', re.MULTILINE)    

    valid     = []
    mismatch  = []

    mycmp     = cmp  # create local cmp for faster calling in the loop

    with open(infile, 'rb') as f:

        alignment = f.read()

    # get sequences and remove header and newlines
    seq_1       = iter(re_1.search(alignment).group(1).replace('\n',''))
    seq_2       = iter(re_2.search(alignment).group(1).replace('\n',''))

    # GENERATOR BLOCK:
    rm_gaps     = ((i, j) for (i, j) in it.izip(seq_1, seq_2) if i != '-')  # remove gaps
    all_sites   = enumerate(rm_gaps, 1)  # enumerate (1-based)
    valid_sites = ((pos, i, j) for (pos, (i, j)) in all_sites if set(i+j).issubset('ACGT'))  # all sites where both sequences are valid

    for (pos, i, j) in valid_sites:

        valid += [pos]
        if mycmp(i,j):
            mismatch += [pos]

    return valid, mismatch

      

        
        
        
      

    

EDIT: per popular request, here is a link to one of the files for people who would like to test their code: https://www.dropbox.com/s/i904fil7cvv1vco/chr1_homo_maca_100Multiz.fa?dl=0