python - 从一组段范围构建掩码

Question

我正在处理由长序列组成的数据(整个人类基因组，因此总序列长度约为 3e9)。我有 22 个文件，每个文件包含 5e7 到 2.5e8 个字符之间的子序列。

就我的问题而言，这些字符要么是0，要么是1，所以文件看起来像这样:

010111101011001001000001100111100011110000110100001011000010

给定的 0 或 1 的索引是它的“位置”(从 0 开始)。

我有一组非重叠范围表示与序列文件中的位置相对应的位置，例如

[(1700, 2000), (9000, 15000), (16000, 18000)]

对于范围列表中的每个位置，我想将序列中的任何 1 转换为 0(转换不包括范围上限，例如一个 python 切片，请参见示例)。

例子:

sequence = 1111011101
ranges   = [(0, 3), (7,10)]
result   = 0001011000  
# the first 3 and last 3 positions are converted to 0 if they are not 
# 0 already, otherwise they are left alone

我正在寻找一种有效的方法来更新给定一组范围的序列，可能多次。我可能会一遍又一遍地这样做，所以我很关心速度。内存不是问题，因此从与序列长度相同的范围创建掩码就可以了，只要创建掩码很快即可。

Answer 1

序列是如何表示的？当然不是这里的整数(前导零呢？)。它是一串数字字符吗？还是整数数组？

假设是整型数组，简单的for循环就没有问题

In [50]: sequence = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0, 1])

In [51]: ranges   = [(0, 3), (7,10)]

In [52]: for r in ranges:
   ....:     sequence[r[0]:r[1]] = 0
   ....:     

In [53]: sequence
Out[53]: array([0, 0, 0, 1, 0, 1, 1, 0, 0, 0])

仅将整个切片设置为零的矢量化和广播操作几乎总是比检查条目是否首先为零的任何操作更快。

如果范围的数量非常大，Python 循环可能会很慢，在这种情况下，您可以简单地将其移动到 Cython，或者考虑并发访问具有 Cython 类型内存 View 的共享内存数组，尤其是如果您可以保证范围永远不会重叠。

如果您从 Python 字符串开始，您可以考虑为它预先计算一个数组格式，例如使用 numpy.char.array。默认情况下，这些数组是不可变的，就像 Python 字符串一样，但您可以将 write 标志设置为 True 以改变它们。如果由于数据序列大小而导致空间成为问题，您可以进一步将数据预先计算为自定义 1 位整数类型的数组，但不要进行此优化，除非某些基准测试表明您确实需要这样做。

假设您可以在 NumPy 中转换为标准长度为 1 的字符串类型，这也有效:

In [69]: s2 = np.char.array("1111011101", itemsize=1)    

In [70]: s2.setflags(write=True)

In [71]: for r in ranges:
    s2[r[0]:r[1]] = '0'
   ....:     

In [72]: s2
Out[72]: 
chararray(['0', '0', '0', '1', '0', '1', '1', '0', '0', '0'], 
      dtype='|S1')

In [73]: s2.tostring()
Out[73]: '0001011000'