如何压缩两个非常大的字符串并返回匹配和不匹配的索引？

Question

我有一组文本文件，其中包含两个长度相同的非常大的字符集。字符集是DNA序列，所以我将它们称为seq_1和{}，它们一起被称为alignment。文件如下所示：

>HEADER1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>HEADER2
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

>HEADER1下序列1中的可能字符是ACGTN-，而{}下序列2中的可能字符是ACGTN-*。在

我想分析序列并返回两个索引列表，我将分别调用valid和{}。在

valid包含两个序列（“对齐”）中的位置在集合中的所有（基于1的）索引ACGT；mismatch包含对齐中的位置在集合ACGT中但彼此不匹配的所有（基于1的）索引。因此mismatch是{}的子集。在

最后一个条件是，在序列1是"-"的位置，我做不递增索引计数器，因为这些在我使用的坐标系中被认为是“间隙”。在

此示例显示了我的预期输出的对齐方式：

^{pr2}$

我希望改进我当前的代码（如下所示），其中包括对序列的regex提取、压缩和枚举非间隙站点到生成器对象中，然后——主要耗时的步骤——在生成器中循环并填充两个列表。我觉得必须有一个基于数组或itertools的解决方案来解决这个问题，它比通过序列及其索引压缩在一起的for循环更有效，我正在寻找建议。在

代码：

^{3}$

编辑：根据流行的需求，这里有一个链接指向某个文件，供那些想测试代码的人使用： https://www.dropbox.com/s/i904fil7cvv1vco/chr1_homo_maca_100Multiz.fa?dl=0

Answer 1

这是我的版本（但你可能开始向我扔石头混淆）请张贴一些较长的测试数据，我想测试它。。。在

seq1='CT-A-CG'
seq2='CNCTA*G'

import numpy as np
def is_not_gap(a): return 0 if (a=='-') else 1
def is_valid(a): return 1 if (a=='A' or a=='C' or a=='G' or a=='T' ) else 0
def is_mm(t): return 0 if t[0]==t[1] else 1

# find the indexes to be retained
retainx=np.where( np.multiply( map(is_not_gap, seq1), map(is_not_gap, seq2) ) )[0].tolist()

# find the valid ones
valid0=np.where(np.multiply( map( is_valid, seq1),map( is_valid, seq2))[retainx])[0].tolist()

# find the mismatches
mm=np.array(map( is_mm, zip( seq1,seq2)))
mismatch0=set(valid0) & set(np.where(mm[retainx])[0])

结果（从零开始的索引）：

^{pr2}$

（如果你愿意，我可以发布更长、更详细的版本）